Index: trunk/LMDZ.TITAN/README =================================================================== --- trunk/LMDZ.TITAN/README (revision 1787) +++ trunk/LMDZ.TITAN/README (revision 1788) @@ -1326,2 +1326,9 @@ Fixed a bug in calchim and physiq_mod on conversion of tendencies to 3D and molar ratio / mass ratio -> we assume same 2D relative tendency on all lat in calchim and adjust longitudinal variations + +== 26/09/2017 == JVO +Get rid of all the old-generic dummy aerosol scheme ( just left scatterers_h for compilation ) as +the new microphysics for Titan will be plugged in +-> even removed sedimentation ( will be done in the microphysical model ) +Removed : aerave.F aeroptproperties.F90 radii_mod.F90 aerosol_mod.F90 aerave_new.F stokes.F90 +vlz_fi.F callsedim.F iniaerosol.F suaer_corrk.F90 newsedim.F Index: trunk/LMDZ.TITAN/libf/phytitan/aerave.F =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/aerave.F (revision 1787) +++ (revision ) @@ -1,182 +1,0 @@ - subroutine aerave( ndata, - & longdata,epdata,omegdata,gdata, - & longref,epref,temp,nir,longir - & ,epir,omegir,gir,qref ) - - - implicit none - -!================================================================== -! -! Purpose -! ------- -! Calculate mean values of aerosol quantities in each band. -! -! Authors -! ------- -! R. Fournier (1996) -! F. Forget (1996) -! -! Called by -! --------- -! suaer_corrk.F90 -! -! Calls -! ----- -! blackl.F -! -!================================================================== - -! -! R.Fournier 02/1996 -! (modif F.Forget 02/1996) -! le spectre est decoupe en "nir" bandes et cette routine calcule -! les donnees radiatives moyenne sur chaque bande : l'optimisation -! est faite pour une temperature au sol "temp" et une epaisseur -! optique de l'atmosphere "epref" a la longueur d'onde "longref" -! -! dans la version actuelle, les ponderations sont independantes de -! l'epaisseur optique : c'est a dire que "omegir", "gir" -! et "epir/epre" sont independants de "epref". -! en effet les ponderations sont choisies pour une solution exacte -! en couche mince et milieu isotherme. -! -! entree -! -! ndata : taille des champs data -! longdata,epdata,omegdata,gdata : proprietes radiative de l'aerosol -! (longdata longueur d'onde en METRES) -! * longref : longueur d'onde a laquelle l'epaisseur optique -! est connue -! * epref : epaisseur optique a longref -! * temp : temperature choisie pour la ponderation (Planck) -! * nir : nombre d'intervals dans la discretisation spectrale -! du GCM -! * longir : longueurs d'onde definissant ces intervals -! -! sortie -! -! * epir : epaisseur optique moyenne pour chaque interval -! * omegir : "scattering albedo" moyen pour chaque interval -! * gir : "assymetry factor" moyen pour chaque interval -! * qref : extinction coefficient at reference wavelength - -!======================================================================= -! output - - REAL longref - REAL epref - REAL temp - INTEGER nir - REAL*8 longir(nir+1) - REAL epir(nir) - REAL omegir(nir) - REAL gir(nir) - -!======================================================================= - - INTEGER iir,nirmx - PARAMETER (nirmx=100) - INTEGER idata,ndata - - DOUBLE PRECISION tmp1 - REAL tmp2,tmp3 - -!======================================================================= -! input - - REAL emit - REAL totalemit(nirmx) - REAL longdata(ndata),epdata(ndata) - & ,omegdata(ndata),gdata(ndata) - REAL qextcorrdata(ndata) - INTEGER ibande,nbande - PARAMETER (nbande=1000) - - REAL long,deltalong - INTEGER ilong - INTEGER i1,i2 - REAL c1,c2 - REAL factep,qextcorr,omeg,g,qref - - long=longref - - -!======================================================================= -! pre-interpolation - ilong=1 - DO idata=2,ndata - IF (long.gt.longdata(idata)) ilong=idata - ENDDO - i1=ilong - i2=ilong+1 - IF (i2.gt.ndata) i2=ndata - IF (long.lt.longdata(1)) i2=1 - IF (i1.eq.i2) THEN - c1=1.E+0 - c2=0.E+0 - ELSE - c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) - c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) - ENDIF - - qref=c1*epdata(i1)+c2*epdata(i2) - factep=qref/epref - DO idata=1,ndata - qextcorrdata(idata)=epdata(idata)/factep - ENDDO -!======================================================================= - - DO iir=1,nir - - deltalong=(longir(iir+1)-longir(iir)) / nbande - totalemit(iir)=0.E+0 - epir(iir)=0.E+0 - omegir(iir)=0.E+0 - gir(iir)=0.E+0 - - DO ibande=1,nbande - - long=longir(iir) + (ibande-0.5E+0) * deltalong - CALL blackl(DBLE(long),DBLE(temp),tmp1) - emit=REAL(tmp1) - -!======================================================================= -! interpolation - ilong=1 - DO idata=2,ndata - IF (long.gt.longdata(idata)) ilong=idata - ENDDO - i1=ilong - i2=ilong+1 - - IF (i2.gt.ndata) i2=ndata - IF (long.lt.longdata(1)) i2=1 - IF (i1.eq.i2) THEN - c1=1.E+0 - c2=0.E+0 - ELSE - c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) - c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) - ENDIF -!======================================================================= - - qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2) - omeg=c1*omegdata(i1)+c2*omegdata(i2) - g=c1*gdata(i1)+c2*gdata(i2) - - totalemit(iir)=totalemit(iir)+deltalong*emit - epir(iir)=epir(iir)+deltalong*emit*qextcorr - omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr - gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g - - ENDDO - - gir(iir)=gir(iir)/omegir(iir) - omegir(iir)=omegir(iir)/epir(iir) - epir(iir)=epir(iir)/totalemit(iir) - - ENDDO - - return - end Index: trunk/LMDZ.TITAN/libf/phytitan/aerave_new.F =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/aerave_new.F (revision 1787) +++ (revision ) @@ -1,223 +1,0 @@ - SUBROUTINE aerave_new ( ndata, - & longdata,epdata,omegdata,gdata, - & longref,epref,temp,nir,longir - & ,epir,omegir,gir,qref,omegaref ) - - - IMPLICIT NONE -c....................................................................... -c -c R.Fournier 02/1996 -c (modif F.Forget 02/1996) -c le spectre est decoupe en "nir" bandes et cette routine calcule -c les donnees radiatives moyenne sur chaque bande : l'optimisation -c est faite pour une temperature au sol "temp" et une epaisseur -c optique de l'atmosphere "epref" a la longueur d'onde "longref" -c -c dans la version actuelle, les ponderations sont independantes de -c l'epaisseur optique : c'est a dire que "omegir", "gir" -c et "epir/epre" sont independants de "epref". -c en effet les ponderations sont choisies pour une solution exacte -c en couche mince et milieu isotherme. -c -c entree -c -c ndata : taille des champs data -c longdata,epdata,omegdata,gdata : proprietes radiative de l'aerosol -c (longdata longueur d'onde en METRES) -c * longref : longueur d'onde a laquelle l'epaisseur optique -c est connue -c * epref : epaisseur optique a longref -c * temp : temperature choisie pour la ponderation (Planck) -c * nir : nombre d'intervals dans la discretisation spectrale -c du GCM -c * longir : longueurs d'onde definissant ces intervals -c -c sortie -c -c * epir : epaisseur optique moyenne pour chaque interval -c * omegir : "scattering albedo" moyen pour chaque interval -c * gir : "assymetry factor" moyen pour chaque interval -c * qref : extinction coefficient at reference wavelength -c * omegaref : single scat. albedo at reference wavelength -c -c....................................................................... -c - REAL longref - REAL epref - REAL temp - INTEGER nir - REAL*8 longir(nir+1) - REAL epir(nir) - REAL omegir(nir) - REAL gir(nir) -c -c....................................................................... -c - INTEGER iir,nirmx - PARAMETER (nirmx=100) - INTEGER idata,ndata -c -c....................................................................... -c - REAL emit - REAL totalemit(nirmx) - REAL longdata(ndata),epdata(ndata) - & ,omegdata(ndata),gdata(ndata) - REAL qextcorrdata(ndata) - INTEGER ibande,nbande - PARAMETER (nbande=1000) - REAL long,deltalong - INTEGER ilong - INTEGER i1,i2 - REAL c1,c2 - REAL factep,qextcorr,omeg,g - REAL qref,omegaref -c -c....................................................................... -c - DOUBLE PRECISION tmp1 - REAL tmp2,tmp3 -c -c - long=longref - - - !if(nir.eq.27)then - !print*,'long',long - !print*,'longdata',longdata - !print*,'epdata',epdata - !print*,'omegdata',omegdata - !print*,'gdata',gdata - !print*,'data looks aok!' - - !print*,'ndata=',ndata - !print*,'longdata',shape(longdata) - !print*,'epdata',shape(epdata) - !print*,'omegdata',shape(omegdata) - !print*,'gdata',shape(gdata) - ! print*,'longref',longref - !print*,'epref',epref - !print*,'temp',temp - !print*,'nir',nir - !print*,'longir',longir - !print*,'epir',epir - !print*,'omegir',gir - !print*,'qref',qref - !print*,'longir=',longir - !stop - !endif - - -c******************************************************** -c interpolation - ilong=1 - DO idata=2,ndata - IF (long.gt.longdata(idata)) ilong=idata - ENDDO - i1=ilong - i2=ilong+1 - IF (i2.gt.ndata) i2=ndata - IF (long.lt.longdata(1)) i2=1 - IF (i1.eq.i2) THEN - c1=1.E+0 - c2=0.E+0 - ELSE - c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) - c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) - ENDIF -c******************************************************** -c - qref=c1*epdata(i1)+c2*epdata(i2) - omegaref=c1*omegdata(i1)+c2*omegdata(i2) - factep=qref/epref - DO idata=1,ndata - qextcorrdata(idata)=epdata(idata)/factep - ENDDO -c -c....................................................................... -c - DO iir=1,nir -c -c....................................................................... -c - deltalong=(longir(iir+1)-longir(iir)) / nbande - totalemit(iir)=0.E+0 - epir(iir)=0.E+0 - omegir(iir)=0.E+0 - gir(iir)=0.E+0 -c -c....................................................................... -c - DO ibande=1,nbande -c -c....................................................................... -c - long=longir(iir) + (ibande-0.5E+0) * deltalong - CALL blackl(DBLE(long),DBLE(temp),tmp1) - emit=REAL(tmp1) -c -c....................................................................... -c -c******************************************************** -c interpolation - ilong=1 - DO idata=2,ndata - IF (long.gt.longdata(idata)) ilong=idata - ENDDO - i1=ilong - i2=ilong+1 - IF (i2.gt.ndata) i2=ndata - IF (long.lt.longdata(1)) i2=1 - IF (i1.eq.i2) THEN - c1=1.E+0 - c2=0.E+0 - ELSE - c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) - c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) - ENDIF -c******************************************************** -c - qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2) - omeg=c1*omegdata(i1)+c2*omegdata(i2) - g=c1*gdata(i1)+c2*gdata(i2) -c -c....................................................................... -c - totalemit(iir)=totalemit(iir)+deltalong*emit - epir(iir)=epir(iir)+deltalong*emit*qextcorr - omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr - gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g -c -c....................................................................... -c - ENDDO -c -c....................................................................... -c - gir(iir)=gir(iir)/omegir(iir) - omegir(iir)=omegir(iir)/epir(iir) - epir(iir)=epir(iir)/totalemit(iir) -c -c....................................................................... -c - ENDDO -c -c...................................................................... -c -c Diagnostic de controle si on moyenne sur tout le spectre vis ou IR : -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -c tmp2=0.E+0 -c DO iir=1,nir -c tmp2=tmp2+totalemit(iir) -c ENDDO -c tmp3=5.67E-8 * temp**4 -c IF (abs((tmp2-tmp3)/tmp3).gt.0.05E+0) THEN -c PRINT *,'!!!! <---> il manque du Planck (voir moyenne.F)' -c PRINT *,'somme des bandes :',tmp2,'--- Planck:',tmp3 -c ENDIF -c -c...................................................................... -c - RETURN - END Index: trunk/LMDZ.TITAN/libf/phytitan/aeropacity.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/aeropacity.F90 (revision 1787) +++ (revision ) @@ -1,184 +1,0 @@ - Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, & - aerosol,reffrad,QREFvis3d,QREFir3d,tau_col) - - use radinc_h, only : L_TAUMAX,naerkind - use aerosol_mod - USE tracer_h, only: noms - use comcstfi_mod, only: g - use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo, & - pres_bottom_strato,pres_top_strato,obs_tau_col_strato - - 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) -! -! 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(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth - REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius - 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 - - real aerosol0 - - INTEGER l,ig,iq,iaer - - LOGICAL,SAVE :: firstcall=.true. -!$OMP THREADPRIVATE(firstcall) - REAL CBRT - EXTERNAL CBRT - - ! for fixed dust profiles - real expfactor, zp - - ! identify tracers - IF (firstcall) THEN - - write(*,*) "Tracers found in aeropacity:" - - 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_back2lay.ne.0) then - print*,'iaero_back2lay= ',iaero_back2lay - endif - - firstcall=.false. - ENDIF ! of IF (firstcall) - - -! --------------------------------------------------------- - - if (noaero) then - aerosol(1:ngrid,1:nlayer,1)=0.0 ! JVO 2017 : Now iaer = 1 is always dummy co2 for noaero case, since we don't use aerosols - endif - -!================================================================== -! Two-layer aerosols (unknown composition) -! S. Guerlet (2013) -!================================================================== - - if (iaero_back2lay .ne.0) then - iaer=iaero_back2lay -! 1. Initialization - aerosol(1:ngrid,1:nlayer,iaer)=0.0 -! 2. Opacity calculation - DO ig=1,ngrid - DO l=1,nlayer-1 - aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) ) - !! 1. below tropospheric layer: no aerosols - IF (pplev(ig,l) .gt. pres_bottom_tropo) THEN - aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer) - !! 2. tropo layer - ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN - aerosol(ig,l,iaer) = obs_tau_col_tropo*aerosol(ig,l,iaer) - !! 3. linear transition - ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN - expfactor=log(obs_tau_col_strato/obs_tau_col_tropo)/log(pres_bottom_strato/pres_top_tropo) - aerosol(ig,l,iaer)= obs_tau_col_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)*aerosol(ig,l,iaer)/1.5 - !! 4. strato layer - ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .gt. pres_top_strato) THEN - aerosol(ig,l,iaer)= obs_tau_col_strato*aerosol(ig,l,iaer) - !! 5. above strato layer: no aerosols - ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN - aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer) - ENDIF - ENDDO - ENDDO - - ! 3. Re-normalize to observed total column - tau_col(:)=0.0 - DO l=1,nlayer - DO ig=1,ngrid - tau_col(ig) = tau_col(ig) & - + aerosol(ig,l,iaer)/(obs_tau_col_tropo+obs_tau_col_strato) - ENDDO - ENDDO - - DO ig=1,ngrid - DO l=1,nlayer-1 - aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig) - ENDDO - ENDDO - - - end if ! if Two-layer aerosols - - -! -------------------------------------------------------------------------- -! 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 - return - end subroutine aeropacity - Index: trunk/LMDZ.TITAN/libf/phytitan/aeroptproperties.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/aeroptproperties.F90 (revision 1787) +++ (revision ) @@ -1,754 +1,0 @@ - SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad, & - QVISsQREF3d,omegaVIS3d,gVIS3d, & - QIRsQREF3d,omegaIR3d,gIR3d, & - QREFvis3d,QREFir3d)!, & -! omegaREFvis3d,omegaREFir3d) - - use radinc_h, only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind - use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir - use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir - use radcommon_h, only: radiustab,nsize - - implicit none - -! ============================================================= -! Aerosol Optical Properties -! -! Description: -! Compute the scattering parameters in each grid -! box, depending on aerosol grain sizes. Log-normal size -! distribution and Gauss-Legendre integration are used. - -! Parameters: -! Don't forget to set the value of varyingnueff below; If -! the effective variance of the distribution for the given -! aerosol is considered homogeneous in the atmosphere, please -! set varyingnueff(iaer) to .false. Resulting computational -! time will be much better. - -! Authors: J.-B. Madeleine, F. Forget, F. Montmessin -! Slightly modified and converted to F90 by R. Wordsworth (2009) -! Varying nueff section removed by R. Wordsworth for simplicity -! ============================================================== - -! Local variables -! --------------- - - - -! ============================================================= - LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false. -! ============================================================= - -! Min. and max radius of the interpolation grid (in METERS) - REAL, PARAMETER :: refftabmin = 2e-8 !2e-8 -! REAL, PARAMETER :: refftabmax = 35e-6 - REAL, PARAMETER :: refftabmax = 1e-3 -! Log of the min and max variance of the interpolation grid - REAL, PARAMETER :: nuefftabmin = -4.6 - REAL, PARAMETER :: nuefftabmax = 0. -! Number of effective radius of the interpolation grid - INTEGER, PARAMETER :: refftabsize = 200 -! Number of effective variances of the interpolation grid -! INTEGER, PARAMETER :: nuefftabsize = 100 - INTEGER, PARAMETER :: nuefftabsize = 1 -! Interpolation grid indices (reff,nueff) - INTEGER :: grid_i,grid_j -! Intermediate variable - REAL :: var_tmp,var3d_tmp(ngrid,nlayer) -! Bilinear interpolation factors - REAL :: kx,ky,k1,k2,k3,k4 -! Size distribution parameters - REAL :: sizedistk1,sizedistk2 -! Pi! - REAL,SAVE :: pi -!$OMP THREADPRIVATE(pi) -! Variables used by the Gauss-Legendre integration: - INTEGER radius_id,gausind - REAL kint - REAL drad - INTEGER, PARAMETER :: ngau = 10 - REAL weightgaus(ngau),radgaus(ngau) - SAVE weightgaus,radgaus -! DATA weightgaus/.2955242247,.2692667193,.2190863625,.1494513491,.0666713443/ -! DATA radgaus/.1488743389,.4333953941,.6794095682,.8650633666,.9739065285/ - DATA radgaus/0.07652652113350,0.22778585114165, & - 0.37370608871528,0.51086700195146, & - 0.63605368072468,0.74633190646476, & - 0.83911697181213,0.91223442826796, & - 0.96397192726078,0.99312859919241/ - - DATA weightgaus/0.15275338723120,0.14917298659407, & - 0.14209610937519,0.13168863843930, & - 0.11819453196154,0.10193011980823, & - 0.08327674160932,0.06267204829828, & - 0.04060142982019,0.01761400714091/ -!$OMP THREADPRIVATE(radgaus,weightgaus) -! Indices - INTEGER :: i,j,k,l,m,iaer,idomain - INTEGER :: ig,lg,chg - -! Local saved variables -! --------------------- - -! Radius axis of the interpolation grid - REAL,SAVE :: refftab(refftabsize) -! Variance axis of the interpolation grid - REAL,SAVE :: nuefftab(nuefftabsize) -! Volume ratio of the grid - REAL,SAVE :: logvratgrid,vratgrid -! Grid used to remember which calculation is done - LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2) = .false. -!$OMP THREADPRIVATE(refftab,nuefftab,logvratgrid,vratgrid,checkgrid) -! Optical properties of the grid (VISIBLE) - REAL,SAVE :: qsqrefVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) - REAL,SAVE :: qextVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) - REAL,SAVE :: qscatVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) - REAL,SAVE :: omegVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) - REAL,SAVE :: gVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) -!$OMP THREADPRIVATE(qsqrefVISgrid,qextVISgrid,qscatVISgrid,omegVISgrid,gVISgrid) -! Optical properties of the grid (INFRARED) - REAL,SAVE :: qsqrefIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) - REAL,SAVE :: qextIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) - REAL,SAVE :: qscatIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) - REAL,SAVE :: omegIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) - REAL,SAVE :: gIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) -!$OMP THREADPRIVATE(qsqrefIRgrid,qextIRgrid,qscatIRgrid,omegIRgrid,gIRgrid) -! Optical properties of the grid (REFERENCE WAVELENGTHS) - REAL,SAVE :: qrefVISgrid(refftabsize,nuefftabsize,naerkind) - REAL,SAVE :: qscatrefVISgrid(refftabsize,nuefftabsize,naerkind) - REAL,SAVE :: qrefIRgrid(refftabsize,nuefftabsize,naerkind) - REAL,SAVE :: qscatrefIRgrid(refftabsize,nuefftabsize,naerkind) - REAL,SAVE :: omegrefVISgrid(refftabsize,nuefftabsize,naerkind) - REAL,SAVE :: omegrefIRgrid(refftabsize,nuefftabsize,naerkind) -!$OMP THREADPRIVATE(qrefVISgrid,qscatrefVISgrid,qrefIRgrid,qscatrefIRgrid,omegrefVISgrid,& - !$OMP omegrefIRgrid) -! Firstcall - LOGICAL,SAVE :: firstcall = .true. -!$OMP THREADPRIVATE(firstcall) -! Variables used by the Gauss-Legendre integration: - REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2) - REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau) - REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau) -!$OMP THREADPRIVATE(normd,dista,distb) - - REAL,SAVE :: radGAUSa(ngau,naerkind,2) - REAL,SAVE :: radGAUSb(ngau,naerkind,2) -!$OMP THREADPRIVATE(radGAUSa,radGAUSb) - - REAL,SAVE :: qsqrefVISa(L_NSPECTV,ngau,naerkind) - REAL,SAVE :: qrefVISa(ngau,naerkind) - REAL,SAVE :: qsqrefVISb(L_NSPECTV,ngau,naerkind) - REAL,SAVE :: qrefVISb(ngau,naerkind) - REAL,SAVE :: omegVISa(L_NSPECTV,ngau,naerkind) - REAL,SAVE :: omegrefVISa(ngau,naerkind) - REAL,SAVE :: omegVISb(L_NSPECTV,ngau,naerkind) - REAL,SAVE :: omegrefVISb(ngau,naerkind) - REAL,SAVE :: gVISa(L_NSPECTV,ngau,naerkind) - REAL,SAVE :: gVISb(L_NSPECTV,ngau,naerkind) -!$OMP THREADPRIVATE(qsqrefVISa,qrefVISa,qsqrefVISb,qrefVISb,omegVISa, & - !$OMP omegrefVISa,omegVISb,omegrefVISb,gVISa,gVISb) - - REAL,SAVE :: qsqrefIRa(L_NSPECTI,ngau,naerkind) - REAL,SAVE :: qrefIRa(ngau,naerkind) - REAL,SAVE :: qsqrefIRb(L_NSPECTI,ngau,naerkind) - REAL,SAVE :: qrefIRb(ngau,naerkind) - REAL,SAVE :: omegIRa(L_NSPECTI,ngau,naerkind) - REAL,SAVE :: omegrefIRa(ngau,naerkind) - REAL,SAVE :: omegIRb(L_NSPECTI,ngau,naerkind) - REAL,SAVE :: omegrefIRb(ngau,naerkind) - REAL,SAVE :: gIRa(L_NSPECTI,ngau,naerkind) - REAL,SAVE :: gIRb(L_NSPECTI,ngau,naerkind) -!$OMP THREADPRIVATE(qsqrefIRa,qrefIRa,qsqrefIRb,qrefIRb,omegIRa,omegrefIRa,& - !$OMP omegIRb,omegrefIRb,gIRa,gIRb) - - REAL :: radiusm - REAL :: radiusr - -! Inputs -! ------ - - INTEGER :: ngrid,nlayer -! Aerosol effective radius used for radiative transfer (meter) - REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) -! Aerosol effective variance used for radiative transfer (n.u.) - REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) - -! Outputs -! ------- - - REAL,INTENT(OUT) :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind) - REAL,INTENT(OUT) :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind) - REAL,INTENT(OUT) :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind) - - REAL,INTENT(OUT) :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind) - REAL,INTENT(OUT) :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind) - REAL,INTENT(OUT) :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind) - - REAL,INTENT(OUT) :: QREFvis3d(ngrid,nlayer,naerkind) - REAL,INTENT(OUT) :: QREFir3d(ngrid,nlayer,naerkind) - - REAL :: omegaREFvis3d(ngrid,nlayer,naerkind) - REAL :: omegaREFir3d(ngrid,nlayer,naerkind) - - DO iaer = 1, naerkind ! Loop on aerosol kind - IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN -!================================================================== -! If there is one single particle size, optical -! properties of the considered aerosol are homogeneous - DO lg = 1, nlayer - DO ig = 1, ngrid - DO chg = 1, L_NSPECTV - QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1) - omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1) - gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1) - ENDDO - DO chg = 1, L_NSPECTI - QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1) - omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1) - gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1) - ENDDO - QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1) - QREFir3d(ig,lg,iaer)=QREFir(iaer,1) - omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1) - omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1) - ENDDO - ENDDO - - - if (firstcall) then - print*,'Optical prop. of the aerosol are homogenous for:' - print*,'iaer = ',iaer - endif - - ELSE ! Varying effective radius and variance - DO idomain = 1, 2 ! Loop on visible or infrared channel -!================================================================== -! 1. Creating the effective radius and variance grid -! -------------------------------------------------- - IF (firstcall) THEN - -! 1.1 Pi! - pi = 2. * asin(1.e0) - -! 1.2 Effective radius - refftab(1) = refftabmin - refftab(refftabsize) = refftabmax - - logvratgrid = log(refftabmax/refftabmin) / float(refftabsize-1)*3. - vratgrid = exp(logvratgrid) - - do i = 2, refftabsize-1 - refftab(i) = refftab(i-1)*vratgrid**(1./3.) - enddo - -! 1.3 Effective variance - if(nuefftabsize.eq.1)then ! addded by RDW - print*,'Warning: no variance range in aeroptproperties' - nuefftab(1)=0.2 - else - do i = 0, nuefftabsize-1 - nuefftab(i+1) = exp( nuefftabmin + i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) ) - enddo - endif - - firstcall = .false. - ENDIF - -! 1.4 Radius middle point and range for Gauss integration - radiusm=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) + radiustab(iaer,idomain,1)) - radiusr=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) - radiustab(iaer,idomain,1)) - -! 1.5 Interpolating data at the Gauss quadrature points: - DO gausind=1,ngau - drad=radiusr*radgaus(gausind) - radGAUSa(gausind,iaer,idomain)=radiusm-drad - - radius_id=minloc(abs(radiustab(iaer,idomain,:) - (radiusm-drad)),DIM=1) - IF ((radiustab(iaer,idomain,radius_id) - (radiusm-drad)).GT.0) THEN - radius_id=radius_id-1 - ENDIF - IF (radius_id.GE.nsize(iaer,idomain)) THEN - radius_id=nsize(iaer,idomain)-1 - kint = 1. - ELSEIF (radius_id.LT.1) THEN - radius_id=1 - kint = 0. - ELSE - kint = ( (radiusm-drad) - & - radiustab(iaer,idomain,radius_id) ) / & - ( radiustab(iaer,idomain,radius_id+1) - & - radiustab(iaer,idomain,radius_id) ) - ENDIF - IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------------- - DO m=1,L_NSPECTV - qsqrefVISa(m,gausind,iaer)= & - (1-kint)*QVISsQREF(m,iaer,radius_id) + & - kint*QVISsQREF(m,iaer,radius_id+1) - omegVISa(m,gausind,iaer)= & - (1-kint)*omegaVIS(m,iaer,radius_id) + & - kint*omegaVIS(m,iaer,radius_id+1) - gVISa(m,gausind,iaer)= & - (1-kint)*gVIS(m,iaer,radius_id) + & - kint*gVIS(m,iaer,radius_id+1) - ENDDO - qrefVISa(gausind,iaer)= & - (1-kint)*QREFvis(iaer,radius_id) + & - kint*QREFvis(iaer,radius_id+1) - omegrefVISa(gausind,iaer)= & - (1-kint)*omegaREFvis(iaer,radius_id) + & - kint*omegaREFvis(iaer,radius_id+1) - ELSE ! INFRARED DOMAIN ---------------------------------- - DO m=1,L_NSPECTI - qsqrefIRa(m,gausind,iaer)= & - (1-kint)*QIRsQREF(m,iaer,radius_id) + & - kint*QIRsQREF(m,iaer,radius_id+1) - omegIRa(m,gausind,iaer)= & - (1-kint)*omegaIR(m,iaer,radius_id) + & - kint*omegaIR(m,iaer,radius_id+1) - gIRa(m,gausind,iaer)= & - (1-kint)*gIR(m,iaer,radius_id) + & - kint*gIR(m,iaer,radius_id+1) - ENDDO - qrefIRa(gausind,iaer)= & - (1-kint)*QREFir(iaer,radius_id) + & - kint*QREFir(iaer,radius_id+1) - omegrefIRa(gausind,iaer)= & - (1-kint)*omegaREFir(iaer,radius_id) + & - kint*omegaREFir(iaer,radius_id+1) - ENDIF - ENDDO - - DO gausind=1,ngau - drad=radiusr*radgaus(gausind) - radGAUSb(gausind,iaer,idomain)=radiusm+drad - - radius_id=minloc(abs(radiustab(iaer,idomain,:) - & - (radiusm+drad)),DIM=1) - IF ((radiustab(iaer,idomain,radius_id) - & - (radiusm+drad)).GT.0) THEN - radius_id=radius_id-1 - ENDIF - IF (radius_id.GE.nsize(iaer,idomain)) THEN - radius_id=nsize(iaer,idomain)-1 - kint = 1. - ELSEIF (radius_id.LT.1) THEN - radius_id=1 - kint = 0. - ELSE - kint = ( (radiusm+drad) - & - radiustab(iaer,idomain,radius_id) ) / & - ( radiustab(iaer,idomain,radius_id+1) - & - radiustab(iaer,idomain,radius_id) ) - ENDIF - IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------------- - DO m=1,L_NSPECTV - qsqrefVISb(m,gausind,iaer)= & - (1-kint)*QVISsQREF(m,iaer,radius_id) + & - kint*QVISsQREF(m,iaer,radius_id+1) - omegVISb(m,gausind,iaer)= & - (1-kint)*omegaVIS(m,iaer,radius_id) + & - kint*omegaVIS(m,iaer,radius_id+1) - gVISb(m,gausind,iaer)= & - (1-kint)*gVIS(m,iaer,radius_id) + & - kint*gVIS(m,iaer,radius_id+1) - ENDDO - qrefVISb(gausind,iaer)= & - (1-kint)*QREFvis(iaer,radius_id) + & - kint*QREFvis(iaer,radius_id+1) - omegrefVISb(gausind,iaer)= & - (1-kint)*omegaREFvis(iaer,radius_id) + & - kint*omegaREFvis(iaer,radius_id+1) - ELSE ! INFRARED DOMAIN ---------------------------------- - DO m=1,L_NSPECTI - qsqrefIRb(m,gausind,iaer)= & - (1-kint)*QIRsQREF(m,iaer,radius_id) + & - kint*QIRsQREF(m,iaer,radius_id+1) - omegIRb(m,gausind,iaer)= & - (1-kint)*omegaIR(m,iaer,radius_id) + & - kint*omegaIR(m,iaer,radius_id+1) - gIRb(m,gausind,iaer)= & - (1-kint)*gIR(m,iaer,radius_id) + & - kint*gIR(m,iaer,radius_id+1) - ENDDO - qrefIRb(gausind,iaer)= & - (1-kint)*QREFir(iaer,radius_id) + & - kint*QREFir(iaer,radius_id+1) - omegrefIRb(gausind,iaer)= & - (1-kint)*omegaREFir(iaer,radius_id) + & - kint*omegaREFir(iaer,radius_id+1) - ENDIF - ENDDO - -!================================================================== -! CONSTANT NUEFF FROM HERE -!================================================================== - -! 2. Compute the scattering parameters using linear -! interpolation over grain sizes and constant nueff -! --------------------------------------------------- - - DO lg = 1,nlayer - DO ig = 1, ngrid -! 2.1 Effective radius index and kx calculation - var_tmp=reffrad(ig,lg,iaer)/refftabmin - var_tmp=log(var_tmp)*3. - var_tmp=var_tmp/logvratgrid+1. - grid_i=floor(var_tmp) - IF (grid_i.GE.refftabsize) THEN -! WRITE(*,*) 'Warning: particle size in grid box #' -! WRITE(*,*) ig,' is too large to be used by the ' -! WRITE(*,*) 'radiative transfer; please extend the ' -! WRITE(*,*) 'interpolation grid to larger grain sizes.' - grid_i=refftabsize-1 - kx = 1. - ELSEIF (grid_i.LT.1) THEN -! WRITE(*,*) 'Warning: particle size in grid box #' -! WRITE(*,*) ig,' is too small to be used by the ' -! WRITE(*,*) 'radiative transfer; please extend the ' -! WRITE(*,*) 'interpolation grid to smaller grain sizes.' - grid_i=1 - kx = 0. - ELSE - kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) / & - ( refftab(grid_i+1)-refftab(grid_i) ) - ENDIF -! 2.3 Integration - DO j=grid_i,grid_i+1 -! 2.3.1 Check if the calculation has been done - IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN -! 2.3.2 Log-normal dist., r_g and sigma_g are defined -! in [hansen_1974], "Light scattering in planetary -! atmospheres", Space Science Reviews 16 527-610. -! Here, sizedistk1=r_g and sizedistk2=sigma_g^2 - sizedistk2 = log(1.+nueffrad(1,1,iaer)) - sizedistk1 = exp(2.5*sizedistk2) - sizedistk1 = refftab(j) / sizedistk1 - - normd(j,1,iaer,idomain) = 1e-30 - DO gausind=1,ngau - drad=radiusr*radgaus(gausind) - dista(j,1,iaer,idomain,gausind) = LOG((radiusm-drad)/sizedistk1) - dista(j,1,iaer,idomain,gausind) = & - EXP(-dista(j,1,iaer,idomain,gausind) * & - dista(j,1,iaer,idomain,gausind) * & - 0.5e0/sizedistk2)/(radiusm-drad) - dista(j,1,iaer,idomain,gausind) = & - dista(j,1,iaer,idomain,gausind) / & - (sqrt(2e0*pi*sizedistk2)) - - distb(j,1,iaer,idomain,gausind) = LOG((radiusm+drad)/sizedistk1) - distb(j,1,iaer,idomain,gausind) = & - EXP(-distb(j,1,iaer,idomain,gausind) * & - distb(j,1,iaer,idomain,gausind) * & - 0.5e0/sizedistk2)/(radiusm+drad) - distb(j,1,iaer,idomain,gausind) = & - distb(j,1,iaer,idomain,gausind) / & - (sqrt(2e0*pi*sizedistk2)) - - normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) + & - weightgaus(gausind) * & - ( & - distb(j,1,iaer,idomain,gausind) * pi * & - radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) + & - dista(j,1,iaer,idomain,gausind) * pi * & - radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) & - ) - ENDDO - IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------- -! 2.3.3.vis Initialization - qsqrefVISgrid(j,1,:,iaer)=0. - qextVISgrid(j,1,:,iaer)=0. - qscatVISgrid(j,1,:,iaer)=0. - omegVISgrid(j,1,:,iaer)=0. - gVISgrid(j,1,:,iaer)=0. - qrefVISgrid(j,1,iaer)=0. - qscatrefVISgrid(j,1,iaer)=0. - omegrefVISgrid(j,1,iaer)=0. - - DO gausind=1,ngau - DO m=1,L_NSPECTV -! Convolution: - qextVISgrid(j,1,m,iaer) = & - qextVISgrid(j,1,m,iaer) + & - weightgaus(gausind) * & - ( & - qsqrefVISb(m,gausind,iaer) * & - qrefVISb(gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - qsqrefVISa(m,gausind,iaer) * & - qrefVISa(gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - qscatVISgrid(j,1,m,iaer) = & - qscatVISgrid(j,1,m,iaer) + & - weightgaus(gausind) * & - ( & - omegVISb(m,gausind,iaer) * & - qsqrefVISb(m,gausind,iaer) * & - qrefVISb(gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - omegVISa(m,gausind,iaer) * & - qsqrefVISa(m,gausind,iaer) * & - qrefVISa(gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - gVISgrid(j,1,m,iaer) = & - gVISgrid(j,1,m,iaer) + & - weightgaus(gausind) * & - ( & - omegVISb(m,gausind,iaer) * & - qsqrefVISb(m,gausind,iaer) * & - qrefVISb(gausind,iaer) * & - gVISb(m,gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - omegVISa(m,gausind,iaer) * & - qsqrefVISa(m,gausind,iaer) * & - qrefVISa(gausind,iaer) * & - gVISa(m,gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - ENDDO - qrefVISgrid(j,1,iaer) = & - qrefVISgrid(j,1,iaer) + & - weightgaus(gausind) * & - ( & - qrefVISb(gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - qrefVISa(gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - qscatrefVISgrid(j,1,iaer) = & - qscatrefVISgrid(j,1,iaer) + & - weightgaus(gausind) * & - ( & - omegrefVISb(gausind,iaer) * & - qrefVISb(gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - omegrefVISa(gausind,iaer) * & - qrefVISa(gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - ENDDO - - qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) / & - normd(j,1,iaer,idomain) - qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / & - normd(j,1,iaer,idomain) - omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / & - qrefVISgrid(j,1,iaer) - DO m=1,L_NSPECTV - qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / & - normd(j,1,iaer,idomain) - qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / & - normd(j,1,iaer,idomain) - gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) / & - qscatVISgrid(j,1,m,iaer) / & - normd(j,1,iaer,idomain) - - qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / & - qrefVISgrid(j,1,iaer) - omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / & - qextVISgrid(j,1,m,iaer) - ENDDO - ELSE ! INFRARED DOMAIN ---------- -! 2.3.3.ir Initialization - qsqrefIRgrid(j,1,:,iaer)=0. - qextIRgrid(j,1,:,iaer)=0. - qscatIRgrid(j,1,:,iaer)=0. - omegIRgrid(j,1,:,iaer)=0. - gIRgrid(j,1,:,iaer)=0. - qrefIRgrid(j,1,iaer)=0. - qscatrefIRgrid(j,1,iaer)=0. - omegrefIRgrid(j,1,iaer)=0. - - DO gausind=1,ngau - DO m=1,L_NSPECTI -! Convolution: - qextIRgrid(j,1,m,iaer) = & - qextIRgrid(j,1,m,iaer) + & - weightgaus(gausind) * & - ( & - qsqrefIRb(m,gausind,iaer) * & - qrefVISb(gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - qsqrefIRa(m,gausind,iaer) * & - qrefVISa(gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - qscatIRgrid(j,1,m,iaer) = & - qscatIRgrid(j,1,m,iaer) + & - weightgaus(gausind) * & - ( & - omegIRb(m,gausind,iaer) * & - qsqrefIRb(m,gausind,iaer) * & - qrefVISb(gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - omegIRa(m,gausind,iaer) * & - qsqrefIRa(m,gausind,iaer) * & - qrefVISa(gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - gIRgrid(j,1,m,iaer) = & - gIRgrid(j,1,m,iaer) + & - weightgaus(gausind) * & - ( & - omegIRb(m,gausind,iaer) * & - qsqrefIRb(m,gausind,iaer) * & - qrefVISb(gausind,iaer) * & - gIRb(m,gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - omegIRa(m,gausind,iaer) * & - qsqrefIRa(m,gausind,iaer) * & - qrefVISa(gausind,iaer) * & - gIRa(m,gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - ENDDO - qrefIRgrid(j,1,iaer) = & - qrefIRgrid(j,1,iaer) + & - weightgaus(gausind) * & - ( & - qrefIRb(gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - qrefIRa(gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - qscatrefIRgrid(j,1,iaer) = & - qscatrefIRgrid(j,1,iaer) + & - weightgaus(gausind) * & - ( & - omegrefIRb(gausind,iaer) * & - qrefIRb(gausind,iaer) * & - pi*radGAUSb(gausind,iaer,idomain) * & - radGAUSb(gausind,iaer,idomain) * & - distb(j,1,iaer,idomain,gausind) + & - omegrefIRa(gausind,iaer) * & - qrefIRa(gausind,iaer) * & - pi*radGAUSa(gausind,iaer,idomain) * & - radGAUSa(gausind,iaer,idomain) * & - dista(j,1,iaer,idomain,gausind) & - ) - ENDDO - - qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) / & - normd(j,1,iaer,idomain) - qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / & - normd(j,1,iaer,idomain) - omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / & - qrefIRgrid(j,1,iaer) - DO m=1,L_NSPECTI - qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / & - normd(j,1,iaer,idomain) - qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / & - normd(j,1,iaer,idomain) - gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) / & - qscatIRgrid(j,1,m,iaer) / & - normd(j,1,iaer,idomain) - - qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / & - qrefVISgrid(j,1,iaer) - omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / & - qextIRgrid(j,1,m,iaer) - ENDDO - ENDIF ! -------------------------- - checkgrid(j,1,iaer,idomain) = .true. - ENDIF !checkgrid - ENDDO !grid_i -! 2.4 Linear interpolation - k1 = (1-kx) - k2 = kx - IF (idomain.EQ.1) THEN ! VISIBLE ------------------------ - DO m=1,L_NSPECTV - QVISsQREF3d(ig,lg,m,iaer) = & - k1*qsqrefVISgrid(grid_i,1,m,iaer) + & - k2*qsqrefVISgrid(grid_i+1,1,m,iaer) - omegaVIS3d(ig,lg,m,iaer) = & - k1*omegVISgrid(grid_i,1,m,iaer) + & - k2*omegVISgrid(grid_i+1,1,m,iaer) - gVIS3d(ig,lg,m,iaer) = & - k1*gVISgrid(grid_i,1,m,iaer) + & - k2*gVISgrid(grid_i+1,1,m,iaer) - ENDDO !L_NSPECTV - QREFvis3d(ig,lg,iaer) = & - k1*qrefVISgrid(grid_i,1,iaer) + & - k2*qrefVISgrid(grid_i+1,1,iaer) - omegaREFvis3d(ig,lg,iaer) = & - k1*omegrefVISgrid(grid_i,1,iaer) + & - k2*omegrefVISgrid(grid_i+1,1,iaer) - ELSE ! INFRARED ----------------------- - DO m=1,L_NSPECTI - QIRsQREF3d(ig,lg,m,iaer) = & - k1*qsqrefIRgrid(grid_i,1,m,iaer) + & - k2*qsqrefIRgrid(grid_i+1,1,m,iaer) - omegaIR3d(ig,lg,m,iaer) = & - k1*omegIRgrid(grid_i,1,m,iaer) + & - k2*omegIRgrid(grid_i+1,1,m,iaer) - gIR3d(ig,lg,m,iaer) = & - k1*gIRgrid(grid_i,1,m,iaer) + & - k2*gIRgrid(grid_i+1,1,m,iaer) - ENDDO !L_NSPECTI - QREFir3d(ig,lg,iaer) = & - k1*qrefIRgrid(grid_i,1,iaer) + & - k2*qrefIRgrid(grid_i+1,1,iaer) - omegaREFir3d(ig,lg,iaer) = & - k1*omegrefIRgrid(grid_i,1,iaer) + & - k2*omegrefIRgrid(grid_i+1,1,iaer) - ENDIF ! -------------------------------- - ENDDO !nlayer - ENDDO !ngrid - -!================================================================== - - - - ENDDO ! idomain - - ENDIF ! nsize = 1 - - ENDDO ! iaer (loop on aerosol kind) - - RETURN - END SUBROUTINE aeroptproperties - - - - Index: trunk/LMDZ.TITAN/libf/phytitan/aerosol_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/aerosol_mod.F90 (revision 1787) +++ (revision ) @@ -1,18 +1,0 @@ -!================================================================== -module aerosol_mod -implicit none -save -!================================================================== - -! aerosol indexes: these are initialized to be 0 if the -! corresponding aerosol was not activated in callphys.def -! -- otherwise a value is given in iniaerosol - logical :: noaero = .false. - -! two-layer simple aerosol model - integer :: iaero_back2lay = 0 -!$OMP THREADPRIVATE(noaero,iaero_back2lay) - -!================================================================== -end module aerosol_mod -!================================================================== Index: trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 (revision 1788) @@ -1,10 +1,10 @@ subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, & albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & - tsurf,fract,dist_star,aerosol, & + tsurf,fract,dist_star, & dtlw,dtsw,fluxsurf_lw, & fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, & fluxabs_sw,fluxtop_dn, & OLR_nu,OSR_nu, & - tau_col,firstcall,lastcall) + firstcall,lastcall) use mod_phys_lmdz_para, only : is_master @@ -14,6 +14,4 @@ use ioipsl_getin_p_mod, only: getin_p use gases_h - use radii_mod, only : su_aer_radii,back2lay_reffrad - use aerosol_mod, only : iaero_back2lay USE tracer_h use comcstfi_mod, only: pi, mugaz, cpp @@ -64,5 +62,4 @@ ! OUTPUT - REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau (kg/kg). REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation. REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation. @@ -75,23 +72,7 @@ REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1). REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1). - REAL,INTENT(OUT) :: tau_col(ngrid) ! Diagnostic from aeropacity. REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015 - ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h. - REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind) - REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind) - REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind) - REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind) - REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind) - REAL :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind) - -! REAL :: omegaREFvis3d(ngrid,nlayer,naerkind) -! REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these... - - REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m) - REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance -!$OMP THREADPRIVATE(reffrad,nueffrad) - !----------------------------------------------------------------------- ! Declaration of the variables required by correlated-k subroutines @@ -114,5 +95,4 @@ REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS) - REAL*8 tauaero(L_LEVELS,naerkind) REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2). @@ -125,5 +105,5 @@ REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo. - INTEGER ig,l,k,nw,iaer + INTEGER ig,l,k,nw real szangle @@ -134,19 +114,6 @@ real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1) - ! Local aerosol optical properties for each column on RADIATIVE grid. - real*8,save,allocatable :: QXVAER(:,:,:) - real*8,save,allocatable :: QSVAER(:,:,:) - real*8,save,allocatable :: GVAER(:,:,:) - real*8,save,allocatable :: QXIAER(:,:,:) - real*8,save,allocatable :: QSIAER(:,:,:) - real*8,save,allocatable :: GIAER(:,:,:) - - real, dimension(:,:,:), save, allocatable :: QREFvis3d - real, dimension(:,:,:), save, allocatable :: QREFir3d -!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER,QREFvis3d,QREFir3d) - ! Miscellaneous : - real*8 temp,temp1,temp2,pweight character(len=10) :: tmp1 character(len=10) :: tmp2 @@ -167,28 +134,6 @@ if(firstcall) then - ! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq) - if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind)) - if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind)) - if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind)) - if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS,L_NSPECTI,naerkind)) - if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS,L_NSPECTI,naerkind)) - if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS,L_NSPECTI,naerkind)) - - !!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...) - IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind)) - IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,naerkind)) - ! Effective radius and variance of the aerosols - IF(.not.ALLOCATED(reffrad)) allocate(reffrad(ngrid,nlayer,naerkind)) - IF(.not.ALLOCATED(nueffrad)) allocate(nueffrad(ngrid,nlayer,naerkind)) - call system('rm -f surf_vals_long.out') - if(naerkind.gt.4)then - print*,'Code not general enough to deal with naerkind > 4 yet.' - call abort - endif - call su_aer_radii(ngrid,nlayer,reffrad,nueffrad) - - !-------------------------------------------------- ! Set up correlated k @@ -207,6 +152,4 @@ call setspv ! Basic visible properties. call sugas_corrk ! Set up gaseous absorption properties. - call suaer_corrk ! Set up aerosol optical properties. - OLR_nu(:,:) = 0. @@ -240,24 +183,4 @@ !======================================================================= - qxvaer(:,:,:)=0.0 - qsvaer(:,:,:)=0.0 - gvaer(:,:,:) =0.0 - - qxiaer(:,:,:)=0.0 - qsiaer(:,:,:)=0.0 - giaer(:,:,:) =0.0 - -!-------------------------------------------------- -! Effective radius and variance of the aerosols -!-------------------------------------------------- - - do iaer=1,naerkind - - if(iaer.eq.iaero_back2lay)then - call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev) - endif - - end do !iaer=1,naerkind. - ! How much light do we get ? @@ -266,16 +189,4 @@ end do - ! Get 3D aerosol optical properties. - call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, & - QVISsQREF3d,omegaVIS3d,gVIS3d, & - QIRsQREF3d,omegaIR3d,gIR3d, & - QREFvis3d,QREFir3d) - - ! Get aerosol optical depths. - call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, & - reffrad,QREFvis3d,QREFir3d, & - tau_col) - - !----------------------------------------------------------------------- @@ -288,139 +199,4 @@ !======================================================================= - -!----------------------------------------------------------------------- -! Aerosol optical properties Qext, Qscat and g. -! The transformation in the vertical is the same as for temperature. -!----------------------------------------------------------------------- - - - do iaer=1,naerkind - ! Shortwave. - do nw=1,L_NSPECTV - - do l=1,nlayer - - temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) & - *QREFvis3d(ig,nlayer+1-l,iaer) - - temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) & - *QREFvis3d(ig,max(nlayer-l,1),iaer) - - qxvaer(2*l,nw,iaer) = temp1 - qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 - - temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer) - temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer) - - qsvaer(2*l,nw,iaer) = temp1 - qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 - - temp1=gvis3d(ig,nlayer+1-l,nw,iaer) - temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer) - - gvaer(2*l,nw,iaer) = temp1 - gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 - - end do ! nlayer - - qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer) - qxvaer(2*nlayer+1,nw,iaer)=0. - - qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer) - qsvaer(2*nlayer+1,nw,iaer)=0. - - gvaer(1,nw,iaer)=gvaer(2,nw,iaer) - gvaer(2*nlayer+1,nw,iaer)=0. - - end do ! L_NSPECTV - - do nw=1,L_NSPECTI - ! Longwave - do l=1,nlayer - - temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) & - *QREFir3d(ig,nlayer+1-l,iaer) - - temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) & - *QREFir3d(ig,max(nlayer-l,1),iaer) - - qxiaer(2*l,nw,iaer) = temp1 - qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 - - temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer) - temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer) - - qsiaer(2*l,nw,iaer) = temp1 - qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 - - temp1=gir3d(ig,nlayer+1-l,nw,iaer) - temp2=gir3d(ig,max(nlayer-l,1),nw,iaer) - - giaer(2*l,nw,iaer) = temp1 - giaer(2*l+1,nw,iaer)=(temp1+temp2)/2 - - end do ! nlayer - - qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer) - qxiaer(2*nlayer+1,nw,iaer)=0. - - qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer) - qsiaer(2*nlayer+1,nw,iaer)=0. - - giaer(1,nw,iaer)=giaer(2,nw,iaer) - giaer(2*nlayer+1,nw,iaer)=0. - - end do ! L_NSPECTI - - end do ! naerkind - - ! Test / Correct for freaky s. s. albedo values. - do iaer=1,naerkind - do k=1,L_LEVELS - - do nw=1,L_NSPECTV - if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then - print*,'Serious problems with qsvaer values' - print*,'in callcorrk' - call abort - endif - if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then - qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer) - endif - end do - - do nw=1,L_NSPECTI - if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then - print*,'Serious problems with qsiaer values' - print*,'in callcorrk' - call abort - endif - if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then - qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer) - endif - end do - - end do ! L_LEVELS - end do ! naerkind - -!----------------------------------------------------------------------- -! Aerosol optical depths -!----------------------------------------------------------------------- - - do iaer=1,naerkind ! a bug was here - do k=0,nlayer-1 - - pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ & - (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1)) - temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer) - tauaero(2*k+2,iaer)=max(temp*pweight,0.d0) - tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0) - - end do - ! boundary conditions - tauaero(1,iaer) = tauaero(2,iaer) - !tauaero(1,iaer) = 0. - - end do ! naerkind ! Albedo and Emissivity. @@ -586,6 +362,5 @@ call optcv(dtauv,tauv,taucumv,plevrad, & - qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, & - tmid,pmid,taugsurf) + wbarv,cosbv,tauray,tmid,pmid,taugsurf) call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, & @@ -629,6 +404,5 @@ call optci(plevrad,tlevrad,dtaui,taucumi, & - qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, & - taugsurfi) + cosbi,wbari,tmid,pmid,taugsurfi) call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, & @@ -655,5 +429,4 @@ print*,'fluxtop_dn=',fluxtop_dn(ig) print*,'acosz=',acosz - print*,'aerosol=',aerosol(ig,:,:) print*,'temp= ',pt(ig,:) print*,'pplay= ',pplay(ig,:) @@ -754,6 +527,4 @@ !$OMP END MASTER !$OMP BARRIER - IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad) - IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad) endif Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 (revision 1788) @@ -31,8 +31,4 @@ logical,save :: mass_redistrib !$OMP THREADPRIVATE(stelbbody,ozone,tracer,mass_redistrib) - logical,save :: sedimentation -!$OMP THREADPRIVATE(sedimentation) - logical,save :: aeroback2lay -!$OMP THREADPRIVATE(aeroback2lay) logical,save :: nosurf logical,save :: oblate @@ -41,24 +37,11 @@ integer,save :: ichim integer,save :: iddist - integer,save :: iaervar integer,save :: iradia integer,save :: startype -!$OMP THREADPRIVATE(ichim,iddist,iaervar,iradia,startype) +!$OMP THREADPRIVATE(ichim,iddist,iradia,startype) real,save :: Fat1AU real,save :: stelTbb !$OMP THREADPRIVATE(Fat1AU,stelTbb) - real,save :: tplanet - real,save :: obs_tau_col_tropo - real,save :: obs_tau_col_strato -!$OMP THREADPRIVATE(tplanet,obs_tau_col_tropo,obs_tau_col_strato) - real,save :: pres_bottom_tropo - real,save :: pres_top_tropo - real,save :: pres_bottom_strato - real,save :: pres_top_strato -!$OMP THREADPRIVATE(pres_bottom_tropo,pres_top_tropo,pres_bottom_strato,pres_top_strato) - real,save :: size_tropo - real,save :: size_strato -!$OMP THREADPRIVATE(size_tropo,size_strato) real,save :: pceil !$OMP THREADPRIVATE(pceil) Index: trunk/LMDZ.TITAN/libf/phytitan/callsedim.F =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/callsedim.F (revision 1787) +++ (revision ) @@ -1,128 +1,0 @@ - SUBROUTINE callsedim(ngrid,nlay, ptimestep, - & pplev,zlev, pt, pdt, - & pq, pdqfi, pdqsed,pdqs_sed,nq) - - use radinc_h, only : naerkind - USE tracer_h, only : radius, rho_q - use comcstfi_mod, only: g - - IMPLICIT NONE - -!================================================================== -! -! Purpose -! ------- -! Calculates sedimentation of aerosols depending on their -! density and radius. -! -! Authors -! ------- -! F. Forget (1999) -! Tracer generalisation by E. Millour (2009) -! -!================================================================== - -c----------------------------------------------------------------------- -c declarations: -c ------------- - -c arguments: -c ---------- - - integer,intent(in):: ngrid ! number of horizontal grid points - integer,intent(in):: nlay ! number of atmospheric layers - real,intent(in):: ptimestep ! physics time step (s) - real,intent(in):: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa) - real,intent(in):: pt(ngrid,nlay) ! temperature at mid-layer (K) - real,intent(in):: pdt(ngrid,nlay) ! tendency on temperature - real,intent(in):: zlev(ngrid,nlay+1) ! altitude at layer boundaries - integer,intent(in) :: nq ! number of tracers - real,intent(in) :: pq(ngrid,nlay,nq) ! tracers (kg/kg) - real,intent(in) :: pdqfi(ngrid,nlay,nq) ! tendency on tracers before - ! sedimentation (kg/kg.s-1) - - real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1) - real,intent(out) :: pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1) - -c local: -c ------ - - INTEGER l,ig, iq - - ! for particles with varying radii: - real reffrad(ngrid,nlay,naerkind) ! particle radius (m) - real nueffrad(ngrid,nlay,naerkind) ! aerosol effective radius variance - - real zqi(ngrid,nlay,nq) ! to locally store tracers - real zt(ngrid,nlay) ! to locally store temperature (K) - real masse (ngrid,nlay) ! Layer mass (kg.m-2) - real epaisseur (ngrid,nlay) ! Layer thickness (m) - real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2) - - - LOGICAL,SAVE :: firstcall=.true. -!$OMP THREADPRIVATE(firstcall) - -c ** un petit test de coherence -c -------------------------- - - IF (firstcall) THEN - firstcall=.false. - ENDIF ! of IF (firstcall) - -!======================================================================= -! Preliminary calculation of the layer characteristics -! (mass (kg.m-2), thickness (m), etc. - - do l=1,nlay - do ig=1, ngrid - masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g - epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l) - zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep - end do - end do - -!====================================================================== -! Calculate the transport due to sedimentation for each tracer -! [This has been rearranged by L. Kerber to allow the sedimentation -! of general tracers.] - - do iq=1,nq - if(radius(iq).gt.1.e-9) then -! (no sedim for gases) - -! store locally updated tracers - - do l=1,nlay - do ig=1, ngrid - zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep - enddo - enddo ! of do l=1,nlay - -!====================================================================== -! Sedimentation -!====================================================================== -! General Case - call newsedim(ngrid,nlay,1,ptimestep, - & pplev,masse,epaisseur,zt,radius(iq),rho_q(iq), - & zqi(1,1,iq),wq) - -!======================================================================= -! Calculate the tendencies -!====================================================================== - - do ig=1,ngrid -! Ehouarn: with new way of tracking tracers by name, this is simply - pdqs_sed(ig,iq) = wq(ig,1)/ptimestep - end do - - DO l = 1, nlay - DO ig=1,ngrid - pdqsed(ig,l,iq)=(zqi(ig,l,iq)- - & (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep - ENDDO - ENDDO - endif ! of no gases - enddo ! of do iq=1,nq - return - end Index: trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90 (revision 1788) @@ -18,7 +18,4 @@ character(len=12),parameter :: surfdir="surface_data" - ! aerdir stores aerosol properties files (optprop_*dat files) - character(LEN=18),parameter :: aerdir="aerosol_properties" - end module datafile_mod !----------------------------------------------------------------------- Index: trunk/LMDZ.TITAN/libf/phytitan/iniaerosol.F =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/iniaerosol.F (revision 1787) +++ (revision ) @@ -1,49 +1,0 @@ - SUBROUTINE iniaerosol() - - - use radinc_h, only: naerkind - use aerosol_mod - use callkeys_mod, only: aeroback2lay - - IMPLICIT NONE -c======================================================================= -c subject: -c -------- -c Initialization related to aerosols -c (Chemical species, ice...) -c -c author: Laura Kerber, S. Guerlet -c ------ -c -c======================================================================= - - integer ia - - ia=0 - - if (aeroback2lay) then - ia=ia+1 - iaero_back2lay=ia - endif - write(*,*) '--- Two-layer aerosol = ', iaero_back2lay - - write(*,*) '=== Number of aerosols= ', ia - -! For the zero aerosol case, we currently make a dummy co2 aerosol which is zero everywhere. -! (See aeropacity.F90 for how this works). A better solution would be to turn off the -! aerosol machinery in the no aerosol case, but this would be complicated. LK - - if (ia.eq.0) then !For the zero aerosol case. - noaero = .true. - ia = 1 - endif - - if (ia.ne.naerkind) then - print*, 'Aerosols counted not equal to naerkind' - print*, 'Compile with tag -s',ia,'to run' - print*, 'or change options in callphys.def' - print*, 'Abort in iniaerosol.F' - call abort - endif - - end Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 (revision 1788) @@ -9,5 +9,5 @@ prad,pg,pr,pcpp) - use radinc_h, only: ini_radinc_h, naerkind + use radinc_h, only: ini_radinc_h use datafile_mod, only: datadir use comdiurn_h, only: sinlat, coslat, sinlon, coslon @@ -385,9 +385,4 @@ endif - write(*,*)"Default planetary temperature?" - tplanet=215.0 - call getin_p("tplanet",tplanet) - write(*,*)" tplanet = ",tplanet - write(*,*)"Which star?" startype=1 ! default value = Sol @@ -399,71 +394,4 @@ call getin_p("Fat1AU",Fat1AU) write(*,*)" Fat1AU = ",Fat1AU - - -! TRACERS: - -! write(*,*)"Number of radiatively active aerosols:" -! naerkind=0. ! default value -! call getin_p("naerkind",naerkind) -! write(*,*)" naerkind = ",naerkind - - -!================================= - - write(*,*)"Radiatively active two-layer aersols?" - aeroback2lay=.false. ! default value - call getin_p("aeroback2lay",aeroback2lay) - write(*,*)" aeroback2lay = ",aeroback2lay - - write(*,*)"TWOLAY AEROSOL: total optical depth ", & - "in the tropospheric layer (visible)" - obs_tau_col_tropo=8.D0 - call getin_p("obs_tau_col_tropo",obs_tau_col_tropo) - write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo - - write(*,*)"TWOLAY AEROSOL: total optical depth ", & - "in the stratospheric layer (visible)" - obs_tau_col_strato=0.08D0 - call getin_p("obs_tau_col_strato",obs_tau_col_strato) - write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato - - write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa" - pres_bottom_tropo=66000.0 - call getin_p("pres_bottom_tropo",pres_bottom_tropo) - write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo - - write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa" - pres_top_tropo=18000.0 - call getin_p("pres_top_tropo",pres_top_tropo) - write(*,*)" pres_top_tropo = ",pres_top_tropo - - write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa" - pres_bottom_strato=2000.0 - call getin_p("pres_bottom_strato",pres_bottom_strato) - write(*,*)" pres_bottom_strato = ",pres_bottom_strato - - write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa" - pres_top_strato=100.0 - call getin_p("pres_top_strato",pres_top_strato) - write(*,*)" pres_top_strato = ",pres_top_strato - - write(*,*)"TWOLAY AEROSOL: particle size in the ", & - "tropospheric layer, in meters" - size_tropo=2.e-6 - call getin_p("size_tropo",size_tropo) - write(*,*)" size_tropo = ",size_tropo - - write(*,*)"TWOLAY AEROSOL: particle size in the ", & - "stratospheric layer, in meters" - size_strato=1.e-7 - call getin_p("size_strato",size_strato) - write(*,*)" size_strato = ",size_strato - -!================================= - - write(*,*) "Gravitationnal sedimentation ?" - sedimentation=.false. ! default value - call getin_p("sedimentation",sedimentation) - write(*,*) " sedimentation = ",sedimentation write(*,*) "Does user want to force cpp and mugaz?" Index: trunk/LMDZ.TITAN/libf/phytitan/initracer.F =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/initracer.F (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/initracer.F (revision 1788) @@ -32,12 +32,7 @@ -c----------------------------------------------------------------------- -c radius(nq) ! aerosol particle radius (m) +c------------------------------------------------ c rho_q(nq) ! tracer densities (kg.m-3) -c qext(nq) ! Single Scat. Extinction coeff at 0.67 um -c alpha_lift(nq) ! saltation vertical flux/horiz flux ratio (m-1) -c alpha_devil(nq) ! lifting coeeficient by dust devil -c----------------------------------------------------------------------- - +c------------------------------------------------ nqtot=nq !! we allocate once for all arrays in common in tracer_h.F90 @@ -45,18 +40,5 @@ IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ALLOCATE(mmol(nq)) - ALLOCATE(radius(nq)) ALLOCATE(rho_q(nq)) - ALLOCATE(qext(nq)) - ALLOCATE(alpha_lift(nq)) - ALLOCATE(alpha_devil(nq)) - ALLOCATE(qextrhor(nq)) - !! initialization - alpha_lift(:)=0. - alpha_devil(:)=0. - - ! Added by JVO 2017 : these arrays are handled later - ! -> initialization is the least we can do, please !!! - radius(:)=0. - qext(:)=0. ! Initialization: get tracer names from the dynamics and check if we are @@ -373,14 +355,3 @@ -c Output for records: -c ~~~~~~~~~~~~~~~~~~ - write(*,*) - Write(*,*) '******** initracer : dust transport parameters :' - write(*,*) 'alpha_lift = ', alpha_lift - write(*,*) 'alpha_devil = ', alpha_devil - write(*,*) 'radius = ', radius - - write(*,*) 'Qext = ', qext - write(*,*) - end Index: trunk/LMDZ.TITAN/libf/phytitan/newsedim.F =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/newsedim.F (revision 1787) +++ (revision ) @@ -1,197 +1,0 @@ - SUBROUTINE newsedim(ngrid,nlay,naersize,ptimestep, - & pplev,masse,epaisseur,pt,rd,rho,pqi,wq) - - use comcstfi_mod, only: r, g - IMPLICIT NONE - -!================================================================== -! -! Purpose -! ------- -! Calculates sedimentation of 1 tracer -! of radius rd (m) and density rho (kg.m-3) -! -!================================================================== - -!----------------------------------------------------------------------- -! declarations -! ------------ - -! arguments -! --------- - - integer,intent(in) :: ngrid ! number of atmospheric columns - integer,intent(in) :: nlay ! number of atmospheric layers - integer,intent(in) :: naersize ! number of particle sizes (1 or number - ! of grid boxes) - real,intent(in) :: ptimestep ! physics time step (s) - real,intent(in) :: pplev(ngrid,nlay+1) ! inter-layer pressures (Pa) - real,intent(in) :: pt(ngrid,nlay) ! mid-layer temperatures (K) - real,intent(in) :: masse (ngrid,nlay) ! atmospheric mass (kg) - real,intent(in) :: epaisseur (ngrid,nlay) ! thickness of atm. layers (m) - real,intent(in) :: rd(naersize) ! particle radius (m) - real,intent(in) :: rho ! particle density (kg.m-3) - real,intent(inout) :: pqi(ngrid,nlay) ! tracer (e.g. ?/kg) - real,intent(out) :: wq(ngrid,nlay+1) ! flux of tracer during timestep (?/m-2) - -c local: -c ------ - - INTEGER l,ig, k, i - REAL rfall - - LOGICAL,SAVE :: firstcall=.true. -!$OMP THREADPRIVATE(firstcall) - -c Traceurs : -c ~~~~~~~~ - real traversee (ngrid,nlay) - real vstokes(ngrid,nlay) - real w(ngrid,nlay) - real ptop, dztop, Ep, Stra - - -c Physical constant -c ~~~~~~~~~~~~~~~~~ -c Gas molecular viscosity (N.s.m-2) - real,parameter :: visc=1.e-5 ! CO2 -c Effective gas molecular radius (m) - real,parameter :: molrad=2.2e-10 ! CO2 - -c local and saved variable - real,save :: a,b -!$OMP THREADPRIVATE(a,b) - -c ** un petit test de coherence -c -------------------------- - - - !print*,'temporary fixed particle rad in newsedim!!' - - IF (firstcall) THEN - firstcall=.false. - - - -!======================================================================= -! Preliminary calculations for sedimenation velocity - -! - Constant to compute stokes speed simple formulae -! (Vstokes = b * rho* r**2 avec b= (2/9) * rho * g / visc - b = 2./9. * g / visc - -! - Constant to compute gas mean free path -! l= (T/P)*a, avec a = ( 0.707*8.31/(4*pi*molrad**2 * avogadro)) - a = 0.707*8.31/(4*3.1416* molrad**2 * 6.023e23) - -c - Correction to account for non-spherical shape (Murphy et al. 1990) -c (correction = 0.85 for irregular particles, 0.5 for disk shaped particles) -c a = a * 0.85 - - - ENDIF - -c----------------------------------------------------------------------- -c 1. initialisation -c ----------------- - -c Sedimentation velocity (m/s) -c ~~~~~~~~~~~~~~~~~~~~~~ -c (stokes law corrected for low pressure by the Cunningham -c slip-flow correction according to Rossow (Icarus 36, 1-50, 1978) - - do l=1,nlay - do ig=1, ngrid - if (naersize.eq.1) then - rfall=rd(1) - else - i=ngrid*(l-1)+ig - rfall=rd(i) ! how can this be correct!!? - endif - - vstokes(ig,l) = b * rho * rfall**2 * - & (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall) - -c Layer crossing time (s) : - traversee(ig,l)= epaisseur(ig,l)/vstokes(ig,l) - end do - end do - - -c Calcul de la masse d'atmosphere correspondant a q transferee -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -c (e.g. on recherche le niveau en dessous de laquelle le traceur -c va traverser le niveau intercouche l : "dztop" est sa hauteur -c au dessus de l (m), "ptop" est sa pression (Pa)) - do l=1,nlay - do ig=1, ngrid - - dztop = vstokes(ig,l)* ptimestep - Ep=0 - k=0 - w(ig,l) = 0. !! JF+AS ajout initialisation (LK MARS) -c ************************************************************** -c Simple Method -cc w(ig,l) = -cc & (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g -cc write(*,*) 'OK simple method l,w =', l, w(ig,l) -cc write(*,*) 'OK simple method dztop =', dztop - w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l))) - !!! Diagnostic: JF. Fix: AS. Date: 05/11 - !!! Probleme arrondi avec la quantite ci-dessus - !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit - !!! ---> dans ce cas on utilise le developpement limite ! - !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2 - - IF ( w(ig,l) .eq. 0. ) THEN - w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g - ELSE - w(ig,l) = w(ig,l) * pplev(ig,l) / g - ENDIF -! LK borrowed simple method from Mars model (AS/JF) - -!************************************************************** -cccc Complex method : - if (dztop.gt.epaisseur(ig,l)) then -cccc Cas ou on "epuise" la couche l : On calcule le flux -cccc Venant de dessus en tenant compte de la variation de Vstokes - - Ep= epaisseur(ig,l) - Stra= traversee(ig,l) - do while(dztop.gt.Ep.and.l+k+1.le.nlay) - k=k+1 - dztop= Ep + vstokes(ig,l+k)*(ptimestep -Stra) - Ep = Ep + epaisseur(ig,l+k) - Stra = Stra + traversee(ig,l+k) - enddo - Ep = Ep - epaisseur(ig,l+k) -! ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) - ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) - IF ( ptop .eq. 1. ) THEN - !PRINT*, 'newsedim: exposant trop petit ', ig, l - ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k))) - ELSE - ptop=pplev(ig,l+k) * ptop - ENDIF - - w(ig,l) = (pplev(ig,l) - ptop)/g - - endif !!! complex method -c -cc write(*,*) 'OK new method l,w =', l, w(ig,l) -cc write(*,*) 'OK new method dztop =', dztop -cc if(l.eq.7)write(*,*)'l=7,k,pplev,Ptop',pplev(ig,l),Ptop -cc if(l.eq.7)write(*,*)'l=7,dztop,Ep',dztop,Ep -cc if(l.eq.6)write(*,*)'l=6,k, w',k, w(1,l) -cc if(l.eq.7)write(*,*)'l=7,k, w',k, w(1,l) -cc if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l) -c ************************************************************** - - end do - end do - - call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq) -c write(*,*) ' newsed: wq(6), wq(7), q(6)', -c & wq(1,6),wq(1,7),pqi(1,6) - - END Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/optci.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/optci.F90 (revision 1788) @@ -1,5 +1,4 @@ subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, & - QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, & - TMID,PMID,TAUGSURF) + COSBI,WBARI,TMID,PMID,TAUGSURF) use radinc_h @@ -41,12 +40,4 @@ real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS) - ! for aerosols - real*8 QXIAER(L_LEVELS,L_NSPECTI,NAERKIND) - real*8 QSIAER(L_LEVELS,L_NSPECTI,NAERKIND) - real*8 GIAER(L_LEVELS,L_NSPECTI,NAERKIND) - real*8 TAUAERO(L_LEVELS,NAERKIND) - real*8 TAUAEROLK(L_LEVELS,L_NSPECTI,NAERKIND) - real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND) - ! Titan customisation ! J. Vatant d'Ollone (2016) @@ -68,5 +59,5 @@ real*8 taugsurf(L_NSPECTI,L_NGAUSS-1) - real*8 DCONT,DAERO + real*8 DCONT double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc double precision p_cross @@ -117,13 +108,4 @@ end do ! levels - ! Spectral dependance of aerosol absorption - do iaer=1,naerkind - DO NW=1,L_NSPECTI - do K=2,L_LEVELS - TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER) - end do ! levels - END DO - end do - do NW=1,L_NSPECTI @@ -132,6 +114,4 @@ ilay = k / 2 ! int. arithmetic => gives the gcm layer index - DAERO=SUM(TAEROS(K,NW,1:naerkind)) ! aerosol absorption - !================= Titan customisation ======================================== call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw)) @@ -240,5 +220,4 @@ DTAUKI(K,nw,ng) = TAUGAS & + DCONT & ! For parameterized continuum absorption - + DAERO & ! For aerosol absorption + DHAZE_T(K,NW) ! For Titan haze @@ -251,5 +230,4 @@ DTAUKI(K,nw,ng) = 0.d0 & + DCONT & ! For parameterized continuum absorption - + DAERO & ! For aerosol absorption + DHAZE_T(K,NW) ! For Titan Haze @@ -262,12 +240,4 @@ ! ====================================================================== - do iaer=1,naerkind - DO NW=1,L_NSPECTI - DO K=2,L_LEVELS - TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) ! effect of scattering albedo - ENDDO - ENDDO - end do - ! Haze scattering DO NW=1,L_NSPECTI @@ -280,6 +250,5 @@ DO L=1,L_NLAYRAD-1 K = 2*L+1 - btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) & - + DHAZES_T(K,NW) + DHAZES_T(K+1,NW) + btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW) END DO ! L vertical loop @@ -287,5 +256,5 @@ L = L_NLAYRAD K = 2*L+1 - btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + DHAZES_T(K,NW) + btemp(L,NW) = DHAZES_T(K,NW) END DO ! NW spectral loop @@ -301,9 +270,4 @@ atemp = 0. if(DTAUI(L,NW,NG) .GT. 1.0D-9) then - do iaer=1,naerkind - atemp = atemp + & - GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + & - GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER) - end do atemp = atemp + & ASF_T(K,NW)*DHAZES_T(K,NW) + & @@ -332,7 +296,4 @@ atemp = 0. if(DTAUI(L,NW,NG) .GT. 1.0D-9) then - do iaer=1,naerkind - atemp = atemp + GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) - end do atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW) WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG) Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 (revision 1788) @@ -1,5 +1,4 @@ SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV, & - QXVAER,QSVAER,GVAER,WBARV,COSBV, & - TAURAY,TAUAERO,TMID,PMID,TAUGSURF) + WBARV,COSBV,TAURAY,TMID,PMID,TAUGSURF) use radinc_h @@ -46,12 +45,4 @@ real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) - - ! for aerosols - real*8 QXVAER(L_LEVELS,L_NSPECTV,NAERKIND) - real*8 QSVAER(L_LEVELS,L_NSPECTV,NAERKIND) - real*8 GVAER(L_LEVELS,L_NSPECTV,NAERKIND) - real*8 TAUAERO(L_LEVELS,NAERKIND) - real*8 TAUAEROLK(L_LEVELS,L_NSPECTV,NAERKIND) - real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND) ! Titan customisation @@ -76,5 +67,5 @@ real*8 taugsurf(L_NSPECTV,L_NGAUSS-1) - real*8 DCONT,DAERO + real*8 DCONT real*8 DRAYAER double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc @@ -129,13 +120,4 @@ end do ! levels - ! Spectral dependance of aerosol absorption - do iaer=1,naerkind - do NW=1,L_NSPECTV - do K=2,L_LEVELS - TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER) - end do ! levels - end do - end do - ! Rayleigh scattering do NW=1,L_NSPECTV @@ -158,8 +140,4 @@ DRAYAER = TRAY(K,NW) ! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity - do iaer=1,naerkind - DRAYAER = DRAYAER + TAEROS(K,NW,IAER) - end do - DRAYAER = DRAYAER + DHAZE_T(K,NW) ! Titan's aerosol @@ -275,12 +253,4 @@ ! we need to calculate the scattering albedo and asymmetry factors - do iaer=1,naerkind - DO NW=1,L_NSPECTV - DO K=2,L_LEVELS - TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) ! effect of scattering albedo - ENDDO - ENDDO - end do - ! Haze scattering DO NW=1,L_NSPECTV @@ -294,8 +264,6 @@ DO L=1,L_NLAYRAD-1 K = 2*L+1 - atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind)) & - + ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW) - btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) & - + DHAZES_T(K,NW) + DHAZES_T(K+1,NW) + atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW) + btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW) ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ? btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW) @@ -306,8 +274,6 @@ L = L_NLAYRAD K = 2*L+1 - atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind)) & - + ASF_T(K,NW)*DHAZES_T(K,NW) - btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) & - + DHAZES_T(K,NW) + atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + btemp(L,NW) = DHAZES_T(K,NW) ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ? btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (revision 1788) @@ -14,5 +14,5 @@ pdu,pdv,pdt,pdq,pdpsrf) - use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind + use radinc_h, only : L_NSPECTI,L_NSPECTV use radcommon_h, only: sigma, glat, grav, BWNV use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe @@ -22,6 +22,5 @@ use geometry_mod, only: latitude, longitude, cell_area USE comgeomfi_h, only: totarea, totarea_planet - USE tracer_h, only: noms, mmol, radius, qext, & - alpha_lift, alpha_devil, qextrhor + USE tracer_h, only: noms, mmol use time_phylmdz_mod, only: ecritphy, iphysiq, nday use phyetat0_mod, only: phyetat0 @@ -74,5 +73,4 @@ ! V. Tracers ! V.1. Chemistry -! V.2. Aerosols and particles. ! V.3. Updates (pressure variations, surface budget). ! V.4. Surface Tracer Update. @@ -225,8 +223,4 @@ ! ----------------- -! Aerosol (dust or ice) extinction optical depth at reference wavelength -! for the "naerkind" optically active aerosols: - - real aerosol(ngrid,nlayer,naerkind) ! Aerosols. real zh(ngrid,nlayer) ! Potential temperature (K). real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!) @@ -274,5 +268,4 @@ real dqsurf(ngrid,nq) ! Cumulated tendencies. real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines. - real zdqssed(ngrid,nq) ! Callsedim routine. real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine. @@ -281,5 +274,4 @@ real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines. real zdqevap(ngrid,nlayer) ! Turbdiff routine. - real zdqsed(ngrid,nlayer,nq) ! Callsedim routine. real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine. @@ -327,11 +319,7 @@ real vmr(ngrid,nlayer) ! volume mixing ratio real time_phys - real,allocatable,dimension(:),save :: tau_col ! Total Aerosol Optical Depth. -!$OMP THREADPRIVATE(tau_col) real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic. - real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2). - ! to test energy conservation (RW+JL) real mass(ngrid,nlayer),massarea(ngrid,nlayer) @@ -353,16 +341,8 @@ real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux. real albedo_equivalent1(ngrid) ! For Equivalent albedo calculation. - real tau_col1(ngrid) ! For aerosol optical depth diagnostic. - real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV) ! For Outgoing Radiation diagnostics. real tf, ntf real,allocatable,dimension(:,:),save :: qsurf_hist !$OMP THREADPRIVATE(qsurf_hist) - - real,allocatable,dimension(:,:,:),save :: reffrad ! Aerosol effective radius (m). By RW - real,allocatable,dimension(:,:,:),save :: nueffrad ! Aerosol effective radius variance. By RW -!$OMP THREADPRIVATE(reffrad,nueffrad) - - real reffcol(ngrid,naerkind) ! Local variables for Titan chemistry and microphysics (JVO 2017) @@ -425,6 +405,4 @@ ALLOCATE(zuprevious(ngrid,nlayer)) ALLOCATE(qsurf_hist(ngrid,nq)) - ALLOCATE(reffrad(ngrid,nlayer,naerkind)) - ALLOCATE(nueffrad(ngrid,nlayer,naerkind)) ALLOCATE(fluxsurf_lw(ngrid)) ALLOCATE(fluxsurf_sw(ngrid)) @@ -439,5 +417,4 @@ ALLOCATE(zdtlw(ngrid,nlayer)) ALLOCATE(zdtsw(ngrid,nlayer)) - ALLOCATE(tau_col(ngrid)) ALLOCATE(dycchi(ngrid,nlayer,nq-nmicro)) ALLOCATE(qysat(nlayer,nq-nmicro)) @@ -455,12 +432,6 @@ dtrad(:,:) = 0.0 fluxrad(:) = 0.0 - tau_col(:) = 0.0 zdtsw(:,:) = 0.0 zdtlw(:,:) = 0.0 - - -! Initialize aerosol indexes. -! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ - call iniaerosol() @@ -781,9 +752,9 @@ call callcorrk(ngrid,nlayer,pq,nq,qsurf, & albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & - tsurf,fract,dist_star,aerosol, & + tsurf,fract,dist_star, & zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, & fluxsurfabs_sw,fluxtop_lw, & fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, & - tau_col,firstcall,lastcall) + firstcall,lastcall) ! Radiative flux from the sky absorbed by the surface (W.m-2). @@ -1057,25 +1028,4 @@ endif - ! ------------------------- - ! V.2. Aerosol particles - ! ------------------------- - - ! Sedimentation. - if (sedimentation) then - - zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0 - zdqssed(1:ngrid,1:nq) = 0.0 - - call callsedim(ngrid,nlayer,ptimestep, & - pplev,zzlev,pt,pdt,pq,pdq, & - zdqsed,zdqssed,nq) - - ! Whether it falls as rain or snow depends only on the surface temperature - pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq) - dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq) - - end if ! end of 'sedimentation' - - ! --------------- ! V.3 Updates @@ -1214,6 +1164,4 @@ print*,'fluxtop_dn has gone crazy' print*,'fluxtop_dn=',fluxtop_dn(ig) - print*,'tau_col=',tau_col(ig) - print*,'aerosol=',aerosol(ig,:,:) print*,'temp= ',pt(ig,:) print*,'pplay= ',pplay(ig,:) @@ -1268,19 +1216,4 @@ call abort endif - - - ! Compute column amounts (kg m-2) if tracers are enabled. - if(tracer)then - qcol(1:ngrid,1:nq)=0.0 - do iq=1,nq - do ig=1,ngrid - qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer)) - enddo - enddo - - ! Generalised for arbitrary aerosols now. By LK - reffcol(1:ngrid,1:naerkind)=0.0 - - endif ! end of 'tracer' @@ -1359,6 +1292,4 @@ call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 'kg m^-2',2,qsurf(1,iq) ) - call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & - 'kg m^-2',2,qcol(1,iq) ) ! call wstats(ngrid,trim(noms(iq))//'_reff', & @@ -1472,18 +1403,7 @@ call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 'kg m^-2',2,qsurf_hist(1,iq) ) - call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & - 'kg m^-2',2,qcol(1,iq) ) ! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & ! 'kg m^-2',2,qsurf(1,iq) ) - - do ig=1,ngrid - if(tau_col(ig).gt.1.e3)then - print*,'WARNING: tau_col=',tau_col(ig) - print*,'at ig=',ig,'in PHYSIQ' - endif - end do - - call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col) enddo ! end of 'nq' loop Index: trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90 (revision 1788) @@ -1,5 +1,4 @@ module radcommon_h - use radinc_h, only: L_NSPECTI, L_NSPECTV, L_NGAUSS, NTstar, NTstop, & - naerkind, nsizemax + use radinc_h, only: L_NSPECTI, L_NSPECTV, L_NGAUSS, NTstar, NTstop implicit none @@ -41,21 +40,4 @@ ! each temperature, pressure, and spectral interval ! -! AEROSOL RADIATIVE OPTICAL CONSTANTS -! -! Shortwave -! ~~~~~~~~~ -! -! For the "naerkind" kind of aerosol radiative properties : -! QVISsQREF : Qext / Qext("longrefvis") -! omegavis : single scattering albedo -! gvis : assymetry factor -! -! Longwave -! ~~~~~~~~ -! -! For the "naerkind" kind of aerosol radiative properties : -! QIRsQREF : Qext / Qext("longrefvis") -! omegaIR : mean single scattering albedo -! gIR : mean assymetry factor REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi @@ -65,8 +47,4 @@ !$OMP WNOV,DWNV,WAVEV,& !$OMP STELLARF,TAURAY) - - REAL*8 blami(L_NSPECTI+1) - REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90 -!$OMP THREADPRIVATE(blami,blamv) !! AS: introduced to avoid doing same computations again for continuum @@ -85,37 +63,4 @@ !$OMP FZEROI,FZEROV) !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk - real QVISsQREF(L_NSPECTV,naerkind,nsizemax) - real omegavis(L_NSPECTV,naerkind,nsizemax) - real gvis(L_NSPECTV,naerkind,nsizemax) - real QIRsQREF(L_NSPECTI,naerkind,nsizemax) - real omegair(L_NSPECTI,naerkind,nsizemax) - real gir(L_NSPECTI,naerkind,nsizemax) -!$OMP THREADPRIVATE(QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir) - - -! Reference wavelengths used to compute reference optical depth (m) -! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - REAL lamrefir(naerkind),lamrefvis(naerkind) - -! Actual number of grain size classes in each domain for a -! given aerosol: - - INTEGER :: nsize(naerkind,2) - -! Particle size axis (depend on the kind of aerosol and the -! radiation domain) - - DOUBLE PRECISION :: radiustab(naerkind,2,nsizemax) -!$OMP THREADPRIVATE(lamrefir,lamrefvis,radiustab) !nsize read by suaer_corrk - -! Extinction coefficient at reference wavelengths; -! These wavelengths are defined in aeroptproperties, and called -! longrefvis and longrefir. - - REAL :: QREFvis(naerkind,nsizemax) - REAL :: QREFir(naerkind,nsizemax) - REAL :: omegaREFvis(naerkind,nsizemax) - REAL :: omegaREFir(naerkind,nsizemax) REAL,SAVE :: tstellar ! Stellar brightness temperature (SW) Index: trunk/LMDZ.TITAN/libf/phytitan/radii_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/radii_mod.F90 (revision 1787) +++ (revision ) @@ -1,121 +1,0 @@ -!================================================================== -module radii_mod -!================================================================== -! module to centralize the radii calculations for aerosols -!================================================================== - - use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo, & - size_tropo,pres_bottom_strato,size_strato - -contains - - -!================================================================== - subroutine su_aer_radii(ngrid,nlayer,reffrad,nueffrad) -!================================================================== -! Purpose -! ------- -! Compute the effective radii of liquid and icy water particles -! -! Authors -! ------- -! Jeremy Leconte (2012) -! -!================================================================== - use ioipsl_getin_p_mod, only: getin_p - use radinc_h, only: naerkind - use aerosol_mod, only: iaero_back2lay - Implicit none - - integer,intent(in) :: ngrid - integer,intent(in) :: nlayer - - real, intent(out) :: reffrad(ngrid,nlayer,naerkind) !aerosols radii (K) - real, intent(out) :: nueffrad(ngrid,nlayer,naerkind) !variance - - logical, save :: firstcall=.true. -!$OMP THREADPRIVATE(firstcall) - integer :: iaer - - print*,'enter su_aer_radii' - do iaer=1,naerkind -! these values will change once the microphysics gets to work -! UNLESS tracer=.false., in which case we should be working with -! a fixed aerosol layer, and be able to define reffrad in a -! .def file. To be improved! - - -! WARNING : Titan adapt. (J. Vatant d'Ollone, 2017) -! - ONLY THE NO AEROSOL CASE FOR NOW SINCE WE COMPUTE THEM ANOTHER WAY ! -! - This routine is just here to keep the code running without unplugging all (yet) -! - There's only the dummy aerosol case on iaer = 1 - if(iaer.eq.1)then - reffrad(1:ngrid,1:nlayer,iaer) = 1.e-4 - nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 - endif - - if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols - reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6 - nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 - endif - - if(iaer.gt.5)then - print*,'Error in callcorrk, naerkind is too high (>5).' - print*,'The code still needs generalisation to arbitrary' - print*,'aerosol kinds and number.' - call abort - endif - - enddo - - print*,'exit su_aer_radii' - - end subroutine su_aer_radii -!================================================================== - - -!================================================================== - subroutine back2lay_reffrad(ngrid,reffrad,nlayer,pplev) -!================================================================== -! Purpose -! ------- -! Compute the effective radii of particles in a 2-layer model -! -! Authors -! ------- -! Sandrine Guerlet (2013) -! -!================================================================== - - use aerosol_mod !! Particle sizes and boundaries of aerosol layers defined there - Implicit none - - integer,intent(in) :: ngrid - - real, intent(out) :: reffrad(ngrid,nlayer) ! particle radii (m) - REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) - INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers - REAL :: expfactor - INTEGER l,ig - - reffrad(:,:)=1e-6 !!initialization, not important - DO ig=1,ngrid - DO l=1,nlayer-1 - IF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN - reffrad(ig,l) = size_tropo - ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN - expfactor=log(size_strato/size_tropo) / log(pres_bottom_strato/pres_top_tropo) - reffrad(ig,l)= size_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor) - ELSEIF (pplev(ig,l) .le. pres_bottom_strato) then - reffrad(ig,l) = size_strato - ENDIF - ENDDO - ENDDO - - end subroutine back2lay_reffrad -!================================================================== - - - -end module radii_mod -!================================================================== Index: trunk/LMDZ.TITAN/libf/phytitan/radinc_h.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/radinc_h.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/radinc_h.F90 (revision 1788) @@ -49,9 +49,4 @@ ! the k-coefficients. Variable component of the mixture ! can in princple be anything: currently it's H2O. -! -! NAERKIND The number of radiatively active aerosol types -! -! NSIZEMAX The maximum number of aerosol particle sizes -! !---------------------------------------------------------------------- @@ -70,5 +65,4 @@ integer, parameter :: L_NSPECTV = NBvisible -! integer, parameter :: NAERKIND = 2 ! set in scatterers.h real, parameter :: L_TAUMAX = 35 @@ -84,9 +78,4 @@ !integer, parameter :: NTstop = 50000 !real*8,parameter :: NTfac = 1.0D+2 - - ! Maximum number of grain size classes for aerosol convolution: - ! This must correspond to size of largest dataset used for aerosol - ! optical properties in datagcm folder. - integer, parameter :: nsizemax = 60 character(len=100),save :: corrkdir Index: trunk/LMDZ.TITAN/libf/phytitan/setspi.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/setspi.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/setspi.F90 (revision 1788) @@ -23,5 +23,5 @@ use radinc_h, only: L_NSPECTI,corrkdir,banddir,NTstar,NTstop,NTfac - use radcommon_h, only: BWNI,BLAMI,WNOI,DWNI,WAVEI,planckir,sigma + use radcommon_h, only: BWNI,WNOI,DWNI,WAVEI,planckir,sigma use datafile_mod, only: datadir use comcstfi_mod, only: pi @@ -136,7 +136,5 @@ DWNI(M) = BWNI(M+1)-BWNI(M) WAVEI(M) = 1.0D+4/WNOI(M) - BLAMI(M) = 0.01D0/BWNI(M) - end do - BLAMI(M) = 0.01D0/BWNI(M) + end do ! note M=L_NSPECTI+1 after loop due to Fortran bizarreness Index: trunk/LMDZ.TITAN/libf/phytitan/setspv.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/setspv.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/setspv.F90 (revision 1788) @@ -24,5 +24,5 @@ use radinc_h, only: L_NSPECTV, corrkdir, banddir - use radcommon_h, only: BWNV,BLAMV,WNOV,DWNV,WAVEV, & + use radcommon_h, only: BWNV,WNOV,DWNV,WAVEV, & STELLARF,TAURAY use datafile_mod, only: datadir @@ -113,7 +113,5 @@ DWNV(M) = BWNV(M+1)-BWNV(M) WAVEV(M) = 1.0E+4/WNOV(M) - BLAMV(M) = 0.01/BWNV(M) end do - BLAMV(M) = 0.01/BWNV(M) ! wavelength in METERS for aerosol stuff ! note M=L_NSPECTV+1 after loop due to Fortran bizarreness Index: trunk/LMDZ.TITAN/libf/phytitan/stokes.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/stokes.F90 (revision 1787) +++ (revision ) @@ -1,68 +1,0 @@ - subroutine stokes(p,t,rd,w,rho_aer) - -!================================================================== -! Purpose -! ------- -! Compute the sedimentation velocity in a low pressure -! atmosphere. -! -! Authors -! ------- -! Francois Forget (1997) -! -!================================================================== - -! use radcommon_h, only : Rgas - use comcstfi_mod, only: g, pi, avocado - - implicit none - -! input -! ----- -! pressure (Pa), Temperature (K), particle radius (m), density - real p, t, rd, rho_aer - -! output -! ------ -! sedimentation velocity (m/s, >0) - real w - -! locally saved variables -! --------------------- - real a,b,molrad,visc - save a,b -!$OMP THREADPRIVATE(a,b) - - LOGICAL firstcall - SAVE firstcall - DATA firstcall/.true./ -!$OMP THREADPRIVATE(firstcall) - - if (firstcall) then - - !print*,'Routine not working: replace Rgas with r' - !stop - !a = 0.707*Rgas/(4*pi*molrad**2 * avocado) - - molrad=2.2e-10 ! CO2 (only used in condense_co2cloud at the moment) - visc=1.e-5 ! CO2 - - - a = 0.707*8.31/(4*pi*molrad**2 * avocado) - b = (2./9.) * rho_aer * g / visc - !print*,'molrad=',molrad - !print*,'visc=',visc - !print*,'a=',a - !print*,'b=',b - !print*,'rho_aer=',rho_aer - !stop - - firstcall=.false. - end if - -! Sedimentation velocity = -! Stokes' Law corrected for low pressures by the Cunningham -! slip-flow correction according to Rossow (Icarus 36, 1-50, 1978) - w = b * rd*rd * (1 + 1.333* (a*T/P)/rd ) - return - end subroutine stokes Index: trunk/LMDZ.TITAN/libf/phytitan/suaer_corrk.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/suaer_corrk.F90 (revision 1787) +++ (revision ) @@ -1,433 +1,0 @@ - subroutine suaer_corrk - - ! inputs - use radinc_h, only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind - use radcommon_h, only: blamv,blami,lamrefir,lamrefvis - use datafile_mod, only: datadir, aerdir - - ! outputs - use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir - use radcommon_h, only: radiustab,nsize,tstellar - use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir - use aerosol_mod - use callkeys_mod, only: tplanet - - implicit none - -!================================================================== -! Purpose -! ------- -! Initialize all radiative aerosol properties -! -! Notes -! ----- -! Reads the optical properties -> Mean -> Variable assignment -! (ASCII files) (see radcommon_h.F90) -! wvl(nwvl) longsun -! ep(nwvl) epav QVISsQREF(L_NSPECTV) -! omeg(nwvl) omegav omegavis(L_NSPECTV) -! gfactor(nwvl) gav gvis(L_NSPECTV) -! -! Authors -! ------- -! Richard Fournier (1996) Francois Forget (1996) -! Frederic Hourdin -! Jean-jacques morcrette *ECMWF* -! MODIF Francois Forget (2000) -! MODIF Franck Montmessin (add water ice) -! MODIF J.-B. Madeleine 2008W27 -! - Optical properties read in ASCII files -! - Add varying radius of the particules -! - Add water ice clouds -! MODIF R. Wordsworth (2009) -! - generalisation to arbitrary spectral bands -! -! WARNING : Titan adapt. (J. Vatant d'Ollone, 2017) -! - ONLY THE NO AEROSOL CASE FOR NOW SINCE WE COMPUTE THEM ANOTHER WAY ! -! - This routine is just here to keep the code running without unplugging all (yet) -!================================================================================================ - -! Optical properties (read in external ASCII files) - INTEGER,SAVE :: nwvl ! Number of wavelengths in - ! the domain (VIS or IR), read by master - -! REAL :: solsir ! visible to infrared ratio -! ! (tau_0.67um/tau_9um) - - REAL, DIMENSION(:),& - ALLOCATABLE, SAVE :: wvl ! Wavelength axis, read by master - REAL, DIMENSION(:),& - ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis, read by master - - REAL, DIMENSION(:,:),& - ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext, read by master - omeg,& ! Single Scattering Albedo, read by master - gfactor ! Assymetry Factor, read by master - -! Local variables: - - INTEGER :: iaer ! Aerosol index - INTEGER :: idomain ! Domain index (1=VIS,2=IR) - INTEGER :: ilw ! longwave index - INTEGER :: isw ! shortwave index - INTEGER :: isize ! Particle size index - INTEGER :: jfile ! ASCII file scan index - INTEGER :: file_unit = 60 - LOGICAL :: file_ok, endwhile - CHARACTER(LEN=132) :: scanline ! ASCII file scanning line - -! I/O of "aerave" (subroutine that spectrally averages -! the single scattering parameters) - - REAL lamref ! reference wavelengths - REAL epref ! reference extinction ep - -! REAL epav(L_NSPECTI) ! Average ep (= /Qext(lamref) if epref=1) -! REAL omegav(L_NSPECTI) ! Average single scattering albedo -! REAL gav(L_NSPECTI) ! Average assymetry parameter - - REAL epavVI(L_NSPECTV) ! Average ep (= /Qext(lamref) if epref=1) - REAL omegavVI(L_NSPECTV) ! Average single scattering albedo - REAL gavVI(L_NSPECTV) ! Average assymetry parameter - - REAL epavIR(L_NSPECTI) ! Average ep (= /Qext(lamref) if epref=1) - REAL omegavIR(L_NSPECTI) ! Average single scattering albedo - REAL gavIR(L_NSPECTI) ! Average assymetry parameter - - logical forgetit ! use francois' data? - integer iwvl - -! Local saved variables: - - CHARACTER(LEN=30), DIMENSION(naerkind,2), SAVE :: file_id -!$OMP THREADPRIVATE(file_id) -!---- Please indicate the names of the optical property files below -! Please also choose the reference wavelengths of each aerosol - - if (noaero) then - print*, 'naerkind= 0' - endif - do iaer=1,naerkind - if (iaer.eq.1) then ! JVO 2017 : Now iaer = 1 is always dummy co2 for noaero case, since we don't use aerosols - forgetit=.true. -! visible - if(forgetit)then - file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values - else - file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat' - endif - lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx ??? - -! infrared - if(forgetit)then - file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values - else - file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat' - endif - lamrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS ??? - endif ! CO2 aerosols -! NOTE: these lamref's are currently for dust, not CO2 ice. -! JB thinks this shouldn't matter too much, but it needs to be -! fixed / decided for the final version. - - if (iaer.eq.iaero_back2lay) then - print*, 'naerkind= back2lay', iaer - -! visible - file_id(iaer,1) = 'optprop_saturn_vis_n20.dat' - lamrefvis(iaer)=0.8E-6 ! -! infrared - file_id(iaer,2) = 'optprop_saturn_ir_n20.dat' - lamrefir(iaer)=6.E-6 ! -! added by SG - endif - - - enddo - -!------------------------------------------------------------------ - -! Initializations: - - radiustab(:,:,:) = 0.0 - gvis(:,:,:) = 0.0 - omegavis(:,:,:) = 0.0 - QVISsQREF(:,:,:) = 0.0 - gir(:,:,:) = 0.0 - omegair(:,:,:) = 0.0 - QIRsQREF(:,:,:) = 0.0 - - DO iaer = 1, naerkind ! Loop on aerosol kind - DO idomain = 1, 2 ! Loop on radiation domain (VIS or IR) -!================================================================== -! 1. READ OPTICAL PROPERTIES -!================================================================== - -! 1.1 Open the ASCII file - -!$OMP MASTER - INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//& - '/'//TRIM(file_id(iaer,idomain)),& - EXIST=file_ok) - IF (file_ok) THEN - OPEN(UNIT=file_unit,& - FILE=TRIM(datadir)//'/'//TRIM(aerdir)//& - '/'//TRIM(file_id(iaer,idomain)),& - FORM='formatted') - ELSE - ! In ye old days these files were stored in datadir; - ! let's be retro-compatible - INQUIRE(FILE=TRIM(datadir)//& - '/'//TRIM(file_id(iaer,idomain)),& - EXIST=file_ok) - IF (file_ok) THEN - OPEN(UNIT=file_unit,& - FILE=TRIM(datadir)//& - '/'//TRIM(file_id(iaer,idomain)),& - FORM='formatted') - ENDIF - ENDIF - IF(.NOT.file_ok) THEN - write(*,*)'suaer_corrk: Problem opening ',& - TRIM(file_id(iaer,idomain)) - write(*,*)'It should be in: ',TRIM(datadir)//'/'//TRIM(aerdir) - write(*,*)'1) You can set this directory address ',& - 'in callphys.def with:' - write(*,*)' datadir = /absolute/path/to/datagcm' - write(*,*)'2) If ',& - TRIM(file_id(iaer,idomain)),& - ' is a LMD reference datafile, it' - write(*,*)' can be obtained online at:' - write(*,*)' http://www.lmd.jussieu.fr/',& - '~lmdz/planets/LMDZ.GENERIC/datagcm/' - CALL ABORT - ENDIF - -! 1.2 Allocate the optical property table - - jfile = 1 - endwhile = .false. - DO WHILE (.NOT.endwhile) - READ(file_unit,*) scanline - IF ((scanline(1:1) .ne. '#').and.& - (scanline(1:1) .ne. ' ')) THEN - BACKSPACE(file_unit) - reading1_seq: SELECT CASE (jfile) ! ==================== - CASE(1) reading1_seq ! nwvl ---------------------------- - read(file_unit,*) nwvl - jfile = jfile+1 - CASE(2) reading1_seq ! nsize --------------------------- - read(file_unit,*) nsize(iaer,idomain) - endwhile = .true. - CASE DEFAULT reading1_seq ! ---------------------------- - WRITE(*,*) 'readoptprop: ',& - 'Error while loading optical properties.' - CALL ABORT - END SELECT reading1_seq ! ============================== - ENDIF - ENDDO - - ALLOCATE(wvl(nwvl)) ! wvl - ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn - ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep - ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg - ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g - - -! 1.3 Read the data - - jfile = 1 - endwhile = .false. - DO WHILE (.NOT.endwhile) - READ(file_unit,*) scanline - IF ((scanline(1:1) .ne. '#').and.& - (scanline(1:1) .ne. ' ')) THEN - BACKSPACE(file_unit) - reading2_seq: SELECT CASE (jfile) ! ==================== - CASE(1) reading2_seq ! wvl ----------------------------- - read(file_unit,*) wvl - jfile = jfile+1 - CASE(2) reading2_seq ! radiusdyn ----------------------- - read(file_unit,*) radiusdyn - jfile = jfile+1 - CASE(3) reading2_seq ! ep ------------------------------ - isize = 1 - DO WHILE (isize .le. nsize(iaer,idomain)) - READ(file_unit,*) scanline - IF ((scanline(1:1) .ne. '#').and.& - (scanline(1:1) .ne. ' ')) THEN - BACKSPACE(file_unit) - read(file_unit,*) ep(:,isize) - isize = isize + 1 - ENDIF - ENDDO - - jfile = jfile+1 - CASE(4) reading2_seq ! omeg ---------------------------- - isize = 1 - DO WHILE (isize .le. nsize(iaer,idomain)) - READ(file_unit,*) scanline - IF ((scanline(1:1) .ne. '#').and.& - (scanline(1:1) .ne. ' ')) THEN - BACKSPACE(file_unit) - read(file_unit,*) omeg(:,isize) - isize = isize + 1 - ENDIF - ENDDO - - jfile = jfile+1 - CASE(5) reading2_seq ! gfactor ------------------------- - isize = 1 - DO WHILE (isize .le. nsize(iaer,idomain)) - READ(file_unit,*) scanline - IF ((scanline(1:1) .ne. '#').and.& - (scanline(1:1) .ne. ' ')) THEN - BACKSPACE(file_unit) - read(file_unit,*) gfactor(:,isize) - isize = isize + 1 - ENDIF - ENDDO - - jfile = jfile+1 - IF ((idomain.NE.1).OR.(iaer.NE.1)) THEN - endwhile = .true. - ENDIF - CASE(6) reading2_seq - endwhile = .true. - CASE DEFAULT reading2_seq ! ---------------------------- - WRITE(*,*) 'readoptprop: ',& - 'Error while loading optical properties.' - CALL ABORT - END SELECT reading2_seq ! ============================== - ENDIF - ENDDO - -! 1.4 Close the file - - CLOSE(file_unit) - -! 1.5 If Francois' data, convert wvl to metres - if(iaer.eq.1)then - if(forgetit)then - ! print*,'please sort out forgetit for naerkind>1' - do iwvl=1,nwvl - wvl(iwvl)=wvl(iwvl)*1.e-6 - enddo - endif - endif - -!$OMP END MASTER -!$OMP BARRIER - - - - - -!================================================================== -! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS -!================================================================== - domain: SELECT CASE (idomain) -!================================================================== - CASE(1) domain ! VISIBLE DOMAIN (idomain=1) -!================================================================== - - lamref=lamrefvis(iaer) - epref=1.E+0 - - DO isize=1,nsize(iaer,idomain) - -! Save the particle sizes - radiustab(iaer,idomain,isize)=radiusdyn(isize) - -! Averaged optical properties (GCM channels) - -! CALL aerave ( nwvl,& -! wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),& -! lamref,epref,tstellar,& -! L_NSPECTV,blamv,epavVI,& -! omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize)) - CALL aerave_new ( nwvl,& - wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),& - lamref,epref,tstellar,& - L_NSPECTV,blamv,epavVI,& - omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize)) - -! Variable assignments (declared in radcommon) - DO isw=1,L_NSPECTV - QVISsQREF(isw,iaer,isize)=epavVI(isw) - gvis(isw,iaer,isize)=gavVI(isw) - omegavis(isw,iaer,isize)=omegavVI(isw) - END DO - - ENDDO - -!================================================================== - CASE(2) domain ! INFRARED DOMAIN (idomain=2) -!================================================================== - - - DO isize=1,nsize(iaer,idomain) ! ---------------------------------- - - lamref=lamrefir(iaer) - epref=1.E+0 - -! Save the particle sizes - radiustab(iaer,idomain,isize)=radiusdyn(isize) - -! Averaged optical properties (GCM channels) - -! epav is /Qext(lamrefir) since epref=1 -! Note: aerave also computes the extinction coefficient at -! the reference wavelength. This is called QREFvis or QREFir -! (not epref, which is a different parameter). -! Reference wavelengths SHOULD BE defined for each aerosol in -! radcommon_h.F90 - -! CALL aerave ( nwvl,& -! wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),& -! lamref,epref,tplanet,& -! L_NSPECTI,blami,epavIR,& -! omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize)) - CALL aerave_new ( nwvl,& - wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),& - lamref,epref,tplanet,& - L_NSPECTI,blami,epavIR,& - omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize)) - - -! Variable assignments (declared in radcommon) - DO ilw=1,L_NSPECTI - QIRsQREF(ilw,iaer,isize)=epavIR(ilw) - gir(ilw,iaer,isize)=gavIR(ilw) - omegair(ilw,iaer,isize)=omegavIR(ilw) - END DO - - - ENDDO ! isize (particle size) ------------------------------------- - - END SELECT domain - -!======================================================================== -! 3. Deallocate temporary variables that were read in the ASCII files -!======================================================================== - -!$OMP BARRIER -!$OMP MASTER - IF (ALLOCATED(wvl)) DEALLOCATE(wvl) ! wvl - IF (ALLOCATED(radiusdyn)) DEALLOCATE(radiusdyn) ! radiusdyn - IF (ALLOCATED(ep)) DEALLOCATE(ep) ! ep - IF (ALLOCATED(omeg)) DEALLOCATE(omeg) ! omeg - IF (ALLOCATED(gfactor)) DEALLOCATE(gfactor) ! g -!$OMP END MASTER -!$OMP BARRIER - - END DO ! Loop on iaer - END DO ! Loop on idomain -!======================================================================== - - RETURN - - - - END subroutine suaer_corrk - Index: trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90 (revision 1787) +++ trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90 (revision 1788) @@ -10,12 +10,7 @@ character*20, DIMENSION(:), ALLOCATABLE :: noms ! name of the tracer real, DIMENSION(:), ALLOCATABLE :: mmol ! mole mass of tracer (g/mol-1) - real, DIMENSION(:), ALLOCATABLE :: radius ! dust and ice particle radius (m) real, DIMENSION(:), ALLOCATABLE :: rho_q ! tracer densities (kg.m-3) - real, DIMENSION(:), ALLOCATABLE :: qext ! Single Scat. Extinction coeff at 0.67 um - real, DIMENSION(:), ALLOCATABLE :: alpha_lift ! saltation vertical flux/horiz flux ratio (m-1) - real, DIMENSION(:), ALLOCATABLE :: alpha_devil ! lifting coeeficient by dust devil - real, DIMENSION(:), ALLOCATABLE :: qextrhor ! Intermediate for computing opt. depth from q -!$OMP THREADPRIVATE(noms,mmol,radius,rho_q,qext,alpha_lift,alpha_devil,qextrhor) +!$OMP THREADPRIVATE(noms,mmol,rho_q) ! tracer indexes: these are initialized in initracer and should be 0 if the Index: trunk/LMDZ.TITAN/libf/phytitan/vlz_fi.F =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/vlz_fi.F (revision 1787) +++ (revision ) @@ -1,190 +1,0 @@ - SUBROUTINE vlz_fi(ngrid,nlayer,q,pente_max,masse,w,wq) -c -c Auteurs: P.Le Van, F.Hourdin, F.Forget -c -c ******************************************************************** -c Shema d'advection " pseudo amont " dans la verticale -c pour appel dans la physique (sedimentation) -c ******************************************************************** -c q rapport de melange (kg/kg)... -c masse : masse de la couche Dp/g -c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2) -c pente_max = 2 conseillee -c -c -c -------------------------------------------------------------------- - IMPLICIT NONE -c -c -c Arguments: -c ---------- - integer,intent(in) :: ngrid, nlayer - real,intent(in) :: masse(ngrid,nlayer),pente_max - REAL,INTENT(INOUT) :: q(ngrid,nlayer) - REAL,INTENT(INOUT) :: w(ngrid,nlayer) - REAL,INTENT(OUT) :: wq(ngrid,nlayer+1) -c -c Local -c --------- -c - INTEGER i,ij,l,j,ii -c - - real dzq(ngrid,nlayer),dzqw(ngrid,nlayer),adzqw(ngrid,nlayer) - real dzqmax - real newmasse - real sigw, Mtot, MQtot - integer m - - REAL SSUM,CVMGP,CVMGT - integer ismax,ismin - - -c On oriente tout dans le sens de la pression c'est a dire dans le -c sens de W - - do l=2,nlayer - do ij=1,ngrid - dzqw(ij,l)=q(ij,l-1)-q(ij,l) - adzqw(ij,l)=abs(dzqw(ij,l)) - enddo - enddo - - do l=2,nlayer-1 - do ij=1,ngrid -#ifdef CRAY - dzq(ij,l)=0.5* - , cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1)) -#else - if(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) then - dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1)) - else - dzq(ij,l)=0. - endif -#endif - dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1)) - dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l)) - enddo - enddo - - do ij=1,ngrid - dzq(ij,1)=0. - dzq(ij,nlayer)=0. - enddo -c --------------------------------------------------------------- -c .... calcul des termes d'advection verticale ....... -c --------------------------------------------------------------- - -c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq -c -c No flux at the model top: - do ij=1,ngrid - wq(ij,nlayer+1)=0. - enddo - -c 1) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION) -c =============================== - - do l = 1,nlayer ! loop different than when w<0 - do ij=1,ngrid - - if(w(ij,l).gt.0.)then - -c Regular scheme (transfered mass < 1 layer) -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if(w(ij,l).le.masse(ij,l))then - sigw=w(ij,l)/masse(ij,l) - wq(ij,l)=w(ij,l)*(q(ij,l)+0.5*(1.-sigw)*dzq(ij,l)) - - -c Extended scheme (transfered mass > 1 layer) -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - else - m=l - Mtot = masse(ij,m) - MQtot = masse(ij,m)*q(ij,m) - if(m.ge.nlayer)goto 88 - do while(w(ij,l).gt.(Mtot+masse(ij,m+1))) - m=m+1 - Mtot = Mtot + masse(ij,m) - MQtot = MQtot + masse(ij,m)*q(ij,m) - if(m.ge.nlayer)goto 88 - end do - 88 continue - if (m.lt.nlayer) then - sigw=(w(ij,l)-Mtot)/masse(ij,m+1) - wq(ij,l)=(MQtot + (w(ij,l)-Mtot)* - & (q(ij,m+1)+0.5*(1.-sigw)*dzq(ij,m+1)) ) - else - w(ij,l) = Mtot - wq(ij,l) = Mqtot - end if - end if - end if - enddo - enddo - -c 2) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION) -c =============================== - goto 99 ! SKIPPING THIS PART FOR SEDIMENTATION - -c Surface flux up: - do ij=1,ngrid - if(w(ij,1).lt.0.) wq(ij,1)=0. ! warning : not always valid - end do - - do l = 1,nlayer-1 ! loop different than when w>0 - do ij=1,ngrid - if(w(ij,l+1).le.0)then - -c Regular scheme (transfered mass < 1 layer) -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if(-w(ij,l+1).le.masse(ij,l))then - sigw=w(ij,l+1)/masse(ij,l) - wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l)) -c Extended scheme (transfered mass > 1 layer) -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - else - m = l-1 - Mtot = masse(ij,m+1) - MQtot = masse(ij,m+1)*q(ij,m+1) - if (m.le.0)goto 77 - do while(-w(ij,l+1).gt.(Mtot+masse(ij,m))) - m=m-1 - Mtot = Mtot + masse(ij,m+1) - MQtot = MQtot + masse(ij,m+1)*q(ij,m+1) - if (m.le.0)goto 77 - end do - 77 continue - - if (m.gt.0) then - sigw=(w(ij,l+1)+Mtot)/masse(ij,m) - wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)* - & (q(ij,m)-0.5*(1.+sigw)*dzq(ij,m)) ) - else -c wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)*qm(ij,1)) - write(*,*) 'a rather weird situation in vlz_fi !' - stop - end if - endif - endif - enddo - enddo - 99 continue - - do l=1,nlayer - do ij=1,ngrid - -cccccccc lines below not used for sedimentation (No real flux) -ccccc newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l) -ccccc q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l)) -ccccc& /newmasse -ccccc masse(ij,l)=newmasse - - q(ij,l)=q(ij,l) + (wq(ij,l+1)-wq(ij,l))/masse(ij,l) - - enddo - enddo - - - end