Index: /trunk/LMDZ.TITAN/libf/phytitan/calc_cpp_mugaz.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/calc_cpp_mugaz.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/calc_cpp_mugaz.F90 (revision 1648) @@ -34,26 +34,15 @@ do igas=1,ngasmx - if(igas.eq.vgas)then - ! ignore variable gas in cpp calculation + ! all values at 300 K from Engineering Toolbox + if(igas.eq.igas_N2)then + mugaz_c = mugaz_c + 28.01*gfrac(igas,nivref) + elseif(igas.eq.igas_H2)then + mugaz_c = mugaz_c + 2.01*gfrac(igas,nivref) + elseif(igas.eq.igas_CH4)then + mugaz_c = mugaz_c + 16.04*gfrac(igas,nivref) else - ! all values at 300 K from Engineering Toolbox - if(igas.eq.igas_N2)then - mugaz_c = mugaz_c + 28.01*gfrac(igas) - elseif(igas.eq.igas_H2)then - mugaz_c = mugaz_c + 2.01*gfrac(igas) - elseif(igas.eq.igas_CH4)then - mugaz_c = mugaz_c + 16.04*gfrac(igas) - elseif(igas.eq.igas_C2H6)then - ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28 - mugaz_c = mugaz_c + 30.07*gfrac(igas) - elseif(igas.eq.igas_C2H2)then - ! C2H2 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1 - mugaz_c = mugaz_c + 26.04*gfrac(igas) - else - print*,'Error in calc_cpp_mugaz: Gas species not recognised!' - call abort - endif + print*,'Error in calc_cpp_mugaz: Gas species not recognised!' + call abort endif - enddo @@ -61,33 +50,16 @@ do igas=1,ngasmx - if(igas.eq.vgas)then - ! ignore variable gas in cpp calculation + ! all values at 300 K from Engineering Toolbox + if(igas.eq.igas_N2)then + cpp_c = cpp_c + 1.040*gfrac(igas,nivref)*28.01/mugaz_c + elseif(igas.eq.igas_H2)then + cpp_c = cpp_c + 14.31*gfrac(igas,nivref)*2.01/mugaz_c + elseif(igas.eq.igas_CH4)then + cpp_c = cpp_c + 2.226*gfrac(igas,nivref)*16.04/mugaz_c else - ! all values at 300 K from Engineering Toolbox - if(igas.eq.igas_N2)then - cpp_c = cpp_c + 1.040*gfrac(igas)*28.01/mugaz_c - elseif(igas.eq.igas_H2)then - cpp_c = cpp_c + 14.31*gfrac(igas)*2.01/mugaz_c - elseif(igas.eq.igas_CH4)then - cpp_c = cpp_c + 2.226*gfrac(igas)*16.04/mugaz_c - elseif(igas.eq.igas_C2H6)then - ! C2H6 - ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28 - cpp_c = cpp_c + 1.763*gfrac(igas)*30.07/mugaz_c - elseif(igas.eq.igas_C2H2)then - ! C2H2 - ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1 - cpp_c = cpp_c + 1.575*gfrac(igas)*26.04/mugaz_c - else - print*,'Error in calc_cpp_mugaz: Gas species not recognised!' - call abort - endif + print*,'Error in calc_cpp_mugaz: Gas species not recognised!' + call abort endif - enddo - - - - cpp_c = 1000.0*cpp_c Index: /trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90 (revision 1648) @@ -26,5 +26,5 @@ use radinc_h, only: L_NSPECTV - use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep + use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, scalep use gases_h use comcstfi_mod, only: g, mugaz @@ -38,5 +38,5 @@ logical typeknown - real*8 tauvar,tauvarvar,tausum,tausumvar,tauwei,tauweivar,bwidth,bstart + real*8 tauvar,tausum,tauwei,bwidth,bstart double precision df @@ -52,8 +52,5 @@ typeknown=.false. do igas=1,ngasmx - if(igas.eq.vgas)then - print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering ' - endif - if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then + if(gfrac(igas,nivref).lt.5.e-2)then print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// & 'as its mixing ratio is less than 0.05.' @@ -74,5 +71,5 @@ endif - if((gfrac(igas).eq.1.0).and.(vgas.eq.0))then + if(gfrac(igas,nivref).eq.1.0)then print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.' typeknown=.true. @@ -92,6 +89,4 @@ tausum = 0.0 tauwei = 0.0 - tausumvar = 0.0 - tauweivar = 0.0 bstart = 10000.0/BWNV(N+1) bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1)) @@ -100,7 +95,6 @@ tauvar=0.0 - tauvarvar=0.0 do igas=1,ngasmx - if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then + if(gfrac(igas,nivref).lt.5.e-2)then ! ignore variable gas in Rayleigh calculation ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation @@ -117,10 +111,5 @@ endif - if(igas.eq.vgas) then - tauvarvar=tauvarvar+tauconsti(igas)*tauvari(igas) - tauvar=tauvar+0.0*0.0*gfrac(igas) - else - tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas) - endif + tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas,nivref) enddo @@ -130,10 +119,7 @@ tauwei=tauwei+df tausum=tausum+tauvar*df - tauweivar=tauweivar+df - tausumvar=tausumvar+tauvarvar*df enddo TAURAY(N)=tausum/tauwei - TAURAYVAR(N)=tausumvar/tauweivar ! we multiply by scalep here (100) because plev, which is used in optcv, ! is in units of mBar, so we need to convert. Index: /trunk/LMDZ.TITAN/libf/phytitan/call_profilgases.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/call_profilgases.F90 (revision 1648) +++ /trunk/LMDZ.TITAN/libf/phytitan/call_profilgases.F90 (revision 1648) @@ -0,0 +1,85 @@ +subroutine call_profilgases(nlayer) + + use gases_h + + implicit none + + + !============================================================================ + ! + ! Purpose + ! ------- + ! Set atmospheric composition for Titan main gases (N2, CH4, H2) from : + ! - Interactive tracers if we run with chemistry + ! or + ! - CH4 vertical profile file (profile.def) at firstcall only + ! + ! Author + ! ------ + ! J. Vatant d'Ollone (2017) + ! + !============================================================================ + + + !-------------------------- + ! 0. Declarations + ! ------------------------- + + integer, intent(in) :: nlayer + + + logical, save :: firstcall=.true. +!$OMP THREADPRIVATE(firstcall) + + + integer ilay, ierr + + + ! ----------------------------------------------- + ! 1. First case : no interactive CH4 tracer + ! ----------------------------------------------- + + if (firstcall) then + +!$OMP MASTER + + ! Load CH4 vertical profile from file 'profile.def' + open(90,file='profile.def',status='old',form='formatted',iostat=ierr) + + if (ierr.eq.0) then + write(*,*) "call_profilgases.F90: reading file profile.def" + read(90,*) ! header + + write(*,*) "call_profilgases.F90: reading CH4 vertical profile..." + + do ilay=1,nlayer + read(90,*,iostat=ierr) gfrac(igas_CH4,ilay) + if (ierr.ne.0) then + write(*,*) 'call_profilgases.F90: error reading CH4 vertical profile... aborting' + call abort + endif + enddo + + ! Then set H2 (fixed) and N2 (what remains) + + gfrac(igas_H2,:) = 1.0E-3 + gfrac(igas_N2,:) = 1.0 - ( gfrac(igas_CH4,:) + gfrac(igas_H2,:) ) + + else + write(*,*) 'Cannot find required file "profile.def"' + call abort + endif + + close(90) +!$OMP END MASTER +!$OMP BARRIER + + firstcall=.false. + endif ! if (firstcall) + + + ! ----------------------------------------------- + ! 2. Second case : interactive CH4 tracer + ! ----------------------------------------------- + +end subroutine call_profilgases Index: /trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 (revision 1648) @@ -1,5 +1,5 @@ subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, & albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & - tsurf,fract,dist_star,aerosol,muvar, & + tsurf,fract,dist_star,aerosol, & dtlw,dtsw,fluxsurf_lw, & fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, & @@ -17,5 +17,5 @@ USE tracer_h use comcstfi_mod, only: pi, mugaz, cpp - use callkeys_mod, only: diurnal,tracer,nosurf,satval, & + use callkeys_mod, only: diurnal,tracer,nosurf, & strictboundcorrk,specOLR @@ -59,5 +59,4 @@ REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day. REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU). - REAL,INTENT(IN) :: muvar(ngrid,nlayer+1) logical,intent(in) :: firstcall ! Signals first call to physics. logical,intent(in) :: lastcall ! Signals last call to physics. @@ -133,5 +132,4 @@ real*8 taugsurf(L_NSPECTV,L_NGAUSS-1) real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1) - real*8 qvar(L_LEVELS) ! Mixing ratio of variable component (mol/mol). ! Local aerosol optical properties for each column on RADIATIVE grid. @@ -153,18 +151,7 @@ character(len=10) :: tmp2 - ! For fixed water vapour profiles. - integer i_var - real RH - real*8 pq_temp(nlayer) -! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!! - real ptemp, Ttemp, qsat - logical OLRz real*8 NFLUXGNDV_nu(L_NSPECTV) - - ! Included by RW for runaway greenhouse 1D study. - real vtmp(nlayer) - REAL*8 muvarrad(L_LEVELS) - + ! Included by MT for albedo calculations. REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation. @@ -459,42 +446,4 @@ else acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude. - endif - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! Note by JL13 : In the following, some indices were changed in the interpolations, -!!! so that the model results are less dependent on the number of layers ! -!!! -!!! --- The older versions are commented with the comment !JL13index --- -!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - do k=1,L_LEVELS - qvar(k) = 1.0D-7 - end do - - - ! IMPORTANT: Now convert from kg/kg to mol/mol. -! do k=1,L_LEVELS -! qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi)) -! end do - - DO l=1,nlayer - muvarrad(2*l) = muvar(ig,nlayer+2-l) - muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2 - END DO - - muvarrad(1) = muvarrad(2) - muvarrad(2*nlayer+1)=muvar(ig,1) - - ! Keep values inside limits for which we have radiative transfer coefficients !!! - if(L_REFVAR.gt.1)then ! (there was a bug here) - do k=1,L_LEVELS - if(qvar(k).lt.wrefvar(1))then - qvar(k)=wrefvar(1)+1.0e-8 - elseif(qvar(k).gt.wrefvar(L_REFVAR))then - qvar(k)=wrefvar(L_REFVAR)-1.0e-8 - endif - end do endif @@ -574,16 +523,16 @@ endif elseif(tlevrad(k).gt.tgasmax)then - print*,'Maximum temperature is outside the radiative' - print*,'transfer kmatrix bounds, exiting.' - print*,"k=",k," tlevrad(k)=",tlevrad(k) - print*,"tgasmax=",tgasmax +! print*,'Maximum temperature is outside the radiative' +! print*,'transfer kmatrix bounds, exiting.' +! print*,"k=",k," tlevrad(k)=",tlevrad(k) +! print*,"tgasmax=",tgasmax if (strictboundcorrk) then call abort else - print*,'***********************************************' - print*,'we allow model to continue with tlevrad=tgasmax' - print*,' ... we assume we know what you are doing ... ' - print*,' ... but do not let this happen too often ... ' - print*,'***********************************************' +! print*,'***********************************************' +! print*,'we allow model to continue with tlevrad=tgasmax' +! print*,' ... we assume we know what you are doing ... ' +! print*,' ... but do not let this happen too often ... ' +! print*,'***********************************************' !tlevrad(k)=tgasmax endif @@ -607,16 +556,16 @@ endif elseif(tmid(k).gt.tgasmax)then - print*,'Maximum temperature is outside the radiative' - print*,'transfer kmatrix bounds, exiting.' - print*,"k=",k," tmid(k)=",tmid(k) - print*,"tgasmax=",tgasmax +! print*,'Maximum temperature is outside the radiative' +! print*,'transfer kmatrix bounds, exiting.' +! print*,"k=",k," tmid(k)=",tmid(k) +! print*,"tgasmax=",tgasmax if (strictboundcorrk) then call abort else - print*,'***********************************************' - print*,'we allow model to continue with tmid=tgasmin' - print*,' ... we assume we know what you are doing ... ' - print*,' ... but do not let this happen too often ... ' - print*,'***********************************************' +! print*,'***********************************************' +! print*,'we allow model to continue with tmid=tgasmin' +! print*,' ... we assume we know what you are doing ... ' +! print*,' ... but do not let this happen too often ... ' +! print*,'***********************************************' tmid(k)=tgasmax endif @@ -649,5 +598,5 @@ call optcv(dtauv,tauv,taucumv,plevrad, & qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, & - tmid,pmid,taugsurf,qvar,muvarrad) + tmid,pmid,taugsurf,gweight) call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, & @@ -692,5 +641,5 @@ call optci(plevrad,tlevrad,dtaui,taucumi, & qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, & - taugsurfi,qvar,muvarrad) + taugsurfi,gweight) call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, & @@ -813,5 +762,4 @@ IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref ) IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref ) - IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar ) IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref ) !$OMP END MASTER Index: /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 (revision 1648) @@ -60,6 +60,5 @@ real,save :: size_tropo real,save :: size_strato - real,save :: satval -!$OMP THREADPRIVATE(size_tropo,size_strato,satval) +!$OMP THREADPRIVATE(size_tropo,size_strato) real,save :: pceil !$OMP THREADPRIVATE(pceil) Index: /trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90 (revision 1648) +++ /trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90 (revision 1648) @@ -0,0 +1,103 @@ +subroutine disr_haze(dz,press,wno,taeros,ssa,cbar) +use datafile_mod, only: datadir +IMPLICIT NONE + +real*8,intent(in) :: dz,press,wno +real*8,intent(inout):: taeros,ssa,cbar + +!--------------------------- +! Attention !! +! taeros is the integrated extinction over the layer +! (extinction * thickness of layer) +!--------------------------- + +integer :: i,j,iw,ip,ierr +real*8 :: wln,factw,factp +integer,parameter :: nbwl_PL=328,nblev_PL=162 +real*8,save :: ext_PL(nblev_PL,nbwl_PL),ssa_PL(nblev_PL,nbwl_PL) +real*8,save :: asf_PL(nblev_PL,nbwl_PL) +real*8,save :: wl_PL(nbwl_PL),press_PL(nblev_PL) +logical,save :: firstcall=.true. +character(len=15) :: dummy + +if (firstcall) then + print*,"DISR HAZE - PANAYOTIS LAVVAS" + + +! read PL table +! wl_PL in nm +! press_PL in Pa + open(11,file=TRIM(datadir)//'/hazetable_PL_original.dat',status="old",iostat=ierr) + read(11,*) dummy,wl_PL + do i=1,nblev_PL + read(11,*) press_PL(i),ext_PL(i,:) ! in cm-1 + enddo + do i=1,nblev_PL + read(11,*) press_PL(i),ssa_PL(i,:) + enddo + do i=1,nblev_PL + read(11,*) press_PL(i),asf_PL(i,:) + enddo + close(11) + ! convert press_PL into millibar for comparison to press in the generic (wtf press in mbar but ok ... to be modified at some point anyway) + press_PL(:)=press_PL(:)*1E-2 + + firstcall=.false. +endif + +! convert wno (in cm-1) into wln (nm) +wln=1E7/wno + +! interpolate the needed values from the table + +iw=1 +do i=2,nbwl_PL + if(wln.gt.wl_PL(i)) then + iw=i + endif +enddo + +ip=1 +do j=2,nblev_PL + if(press.lt.press_PL(j)) then + ip=j + endif +enddo + +if(iw.ne.nbwl_PL) then + factw=(wln-wl_PL(iw))/(wl_PL(iw+1)-wl_PL(iw)) +endif +if(ip.ne.nblev_PL) then + factp=(log10(press)-log10(press_PL(ip)))/(log10(press_PL(ip+1))-log10(press_PL(ip))) +else +! factp=(log10(press)-log10(press_PL(nblev_PL)))/(log10(1.e-10)-log10(press_PL(nblev_PL))) + factp=0. +endif + +if ((iw.ne.nbwl_PL).and.(ip.ne.nblev_PL)) then + taeros= ext_PL(ip,iw) *(1-factw)*(1-factp)+ext_PL(ip+1,iw) *(1-factw)*factp & + +ext_PL(ip,iw+1)*factw *(1-factp)+ext_PL(ip+1,iw+1)*factw *factp + ssa = ssa_PL(ip,iw) *(1-factw)*(1-factp)+ssa_PL(ip+1,iw) *(1-factw)*factp & + +ssa_PL(ip,iw+1)*factw *(1-factp)+ssa_PL(ip+1,iw+1)*factw *factp + cbar = asf_PL(ip,iw) *(1-factw)*(1-factp)+asf_PL(ip+1,iw) *(1-factw)*factp & + +asf_PL(ip,iw+1)*factw *(1-factp)+asf_PL(ip+1,iw+1)*factw *factp +else if ((iw.eq.nbwl_PL).and.(ip.ne.nblev_PL)) then + taeros= ext_PL(ip,iw)*(1-factp)+ext_PL(ip+1,iw)*factp + ssa = ssa_PL(ip,iw)*(1-factp)+ssa_PL(ip+1,iw)*factp + cbar = asf_PL(ip,iw)*(1-factp)+asf_PL(ip+1,iw)*factp +else if ((iw.ne.nbwl_PL).and.(ip.eq.nblev_PL)) then + taeros= ext_PL(ip,iw) *(1-factw)*(1-factp) & + +ext_PL(ip,iw+1)*factw *(1-factp) + ssa = ssa_PL(ip,iw) *(1-factw)*(1-factp) & + +ssa_PL(ip,iw+1)*factw *(1-factp) + cbar = asf_PL(ip,iw) *(1-factw)*(1-factp) & + +asf_PL(ip,iw+1)*factw *(1-factp) +else if ((iw.eq.nbwl_PL).and.(ip.eq.nblev_PL)) then + taeros= ext_PL(ip,iw)*(1-factp) + ssa = ssa_PL(ip,iw)*(1-factp) + cbar = asf_PL(ip,iw)*(1-factp) +endif + +taeros=taeros*dz*1.E2 ! ext in cm-1 * thickness in m * 1E2 + +end subroutine disr_haze Index: /trunk/LMDZ.TITAN/libf/phytitan/gases_h.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/gases_h.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/gases_h.F90 (revision 1648) @@ -3,26 +3,29 @@ implicit none -!======================================================================C +!============================================================================================C ! ! gases_h ! -! A F90-allocatable version for gases.h -- AS 12/2011 +! * A F90-allocatable version for gases.h -- AS 12/2011 ! -!======================================================================C +! * Titan's version : J. Vatant d'Ollone (2017) +! * gfrac is now 2D (altitude-dependent for the CIA-calculations) +! * Reference level (nivref) is added for the cpp_mu and rayleigh routines +!===========================================================================================C ! Set and allocated in su_gases.F90 integer :: ngasmx - integer :: vgas - character*20,allocatable,DIMENSION(:) :: gnom ! name of the gas, read by master - real,allocatable,DIMENSION(:) :: gfrac + integer :: nivref + character*20,allocatable,DIMENSION(:) :: gnom ! name of the gas, set by master + real,allocatable,DIMENSION(:,:) :: gfrac ! in analogy with tracer.h ... - integer :: igas_H2 + integer :: igas_N2 integer :: igas_CH4 - integer :: igas_C2H2 - integer :: igas_C2H6 -!!$OMP THREADPRIVATE(ngasmx,vgas,gnom,gfrac,& -! !$OMP igas_H2,igas_N2,igas_CH4,igas_C2H2,igas_C2H6) + integer :: igas_H2 + +!!$OMP THREADPRIVATE(ngasmx,nivref,gnom,gfrac,& +! !$OMP igas_H2,igas_N2,igas_CH4) end module gases_h Index: /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 (revision 1648) @@ -439,9 +439,4 @@ !================================= - - write(*,*)"What is the saturation % of the variable species?" - satval=0.8 - call getin_p("satval",satval) - write(*,*)" satval = ",satval write(*,*) "Gravitationnal sedimentation ?" @@ -477,7 +472,10 @@ ! mugaz=8.314*1000./pr endif ! of if (force_cpp) - call su_gases + + + call su_gases(nlayer,tracer) call calc_cpp_mugaz - + + PRINT*,'--------------------------------------------' PRINT* @@ -519,5 +517,4 @@ ENDDO - ! initialize variables in radinc_h call ini_radinc_h(nlayer) Index: /trunk/LMDZ.TITAN/libf/phytitan/interpolateCH4CH4.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/interpolateCH4CH4.F90 (revision 1648) +++ /trunk/LMDZ.TITAN/libf/phytitan/interpolateCH4CH4.F90 (revision 1648) @@ -0,0 +1,128 @@ +subroutine interpolateCH4CH4(wn,temp,pres,abcoef,firstcall,ind) + + !================================================================== + ! + ! Purpose + ! ------- + ! Calculates the CH4-CH4 CIA opacity, using a lookup table from + ! HITRAN (2011) + ! + ! Authors + ! ------- + ! J. Vatant (2016) based on R. Wordsworth (2011) + ! + !================================================================== + + use datafile_mod, only: datadir + implicit none + + ! input + double precision wn ! wavenumber (cm^-1) + double precision temp ! temperature (Kelvin) + double precision pres ! pressure (Pascals) + integer :: ind + + ! output + double precision abcoef ! absorption coefficient (m^-1) + + integer nS,nT + parameter(nS=1018) + parameter(nT=10) + + double precision, parameter :: losch = 2.6867774e19 + ! Loschmit's number (molecule cm^-3 at STP) + ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2 + ! see Richard et al. 2011, JQSRT for details + + double precision amagat + double precision wn_arr(nS) + double precision temp_arr(nT) + double precision abs_arr(nS,nT) + + integer k,iT + logical firstcall + + save wn_arr, temp_arr, abs_arr !read by master + + character*100 dt_file + integer strlen,ios + + character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" + + character*20 bleh + double precision blah, Ttemp + integer nres + + + if(temp.gt.400)then + print*,'Your temperatures are too high for this CH4-CH4 CIA dataset.' + print*,'Currently, HITRAN provides data for this pair in the range 40-400 K.' + stop + endif + + amagat = (273.15/temp)*(pres/101325.0) + + if(firstcall)then ! called by sugas_corrk only + print*,'----------------------------------------------------' + print*,'Initialising CH4-CH4 continuum from HITRAN database...' + + ! 1.1 Open the ASCII files + dt_file=TRIM(datadir)//'/continuum_data/CH4-CH4_2011.cia' + +!$OMP MASTER + open(33,file=dt_file,form='formatted',status='old',iostat=ios) + if (ios.ne.0) then ! file not found + write(*,*) 'Error from interpolateCH4CH4' + write(*,*) 'Data file ',trim(dt_file),' not found.' + write(*,*) 'Check that your path to datagcm:',trim(datadir) + write(*,*) 'is correct. You can change it in callphys.def with:' + write(*,*) 'datadir = /absolute/path/to/datagcm' + write(*,*) 'Also check that the continuum data continuum_data/CH4-CH4_2011.cia is there.' + call abort + else + + do iT=1,nT + + read(33,fmat1) bleh,blah,blah,nres,Ttemp + if(nS.ne.nres)then + print*,'Resolution given in file: ',trim(dt_file) + print*,'is ',nres,', which does not match nS.' + print*,'Please adjust nS value in interpolateCH4CH4.F90' + stop + endif + temp_arr(iT)=Ttemp + + do k=1,nS + read(33,*) wn_arr(k),abs_arr(k,it) + end do + + end do + + endif + close(33) +!$OMP END MASTER +!$OMP BARRIER + + print*,'interpolateCH4CH4: At wavenumber ',wn,' cm^-1' + print*,' temperature ',temp,' K' + print*,' pressure ',pres,' Pa' + + endif + call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) + +! print*,'the absorption is ',abcoef,' cm^5 molecule^-2' +! print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' + + abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1 +! abcoef=0. + +! print*,'We have ',amagat,' amagats of CH4' +! print*,'So the absorption is ',abcoef,' m^-1' + +! unlike for Rayleigh scattering, we do not currently weight by the BB function +! however our bands are normally thin, so this is no big deal. + + + return +end subroutine interpolateCH4CH4 + Index: /trunk/LMDZ.TITAN/libf/phytitan/interpolateN2CH4.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/interpolateN2CH4.F90 (revision 1648) +++ /trunk/LMDZ.TITAN/libf/phytitan/interpolateN2CH4.F90 (revision 1648) @@ -0,0 +1,133 @@ +subroutine interpolateN2CH4(wn,temp,presN2,presCH4,abcoef,firstcall,ind) + + !================================================================== + ! + ! Purpose + ! ------- + ! Calculates the N2-CH4 CIA opacity, using a lookup table from + ! HITRAN (2011) + ! + ! Authors + ! ------- + ! J. Vatant (2016) based on R. Wordsworth (2011) + ! + !================================================================== + + use datafile_mod, only: datadir + implicit none + + ! input + double precision wn ! wavenumber (cm^-1) + double precision temp ! temperature (Kelvin) + double precision presN2 ! N2 partial pressure (Pascals) + double precision presCH4 ! CH4 partial pressure (Pascals) + + ! output + double precision abcoef ! absorption coefficient (m^-1) + + integer nS,nT + parameter(nS=1407) + parameter(nT=10) + + double precision, parameter :: losch = 2.6867774e19 + ! Loschmit's number (molecule cm^-3 at STP) + ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2 + ! see Richard et al. 2011, JQSRT for details + + double precision amagatN2 + double precision amagatCH4 + double precision wn_arr(nS) + double precision temp_arr(nT) + double precision abs_arr(nS,nT) + + integer k,iT + logical firstcall + + save wn_arr, temp_arr, abs_arr !read by master + + character*100 dt_file + integer strlen,ios + + character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" + + character*20 bleh + double precision blah, Ttemp + integer nres + integer ind + + if(temp.gt.400)then + print*,'Your temperatures are too high for this N2-CH4 CIA dataset.' + print*,'Please run mixed N2-CH4 atmospheres below T = 400 K.' + stop + endif + + amagatN2 = (273.15/temp)*(presN2/101325.0) + amagatCH4 = (273.15/temp)*(presCH4/101325.0) + + if(firstcall)then ! called by sugas_corrk only + print*,'----------------------------------------------------' + print*,'Initialising N2-CH4 continuum from HITRAN database...' + + ! 1.1 Open the ASCII files + dt_file=TRIM(datadir)//'/continuum_data/N2-CH4_2011.cia' + +!$OMP MASTER + open(33,file=dt_file,form='formatted',status='old',iostat=ios) + if (ios.ne.0) then ! file not found + write(*,*) 'Error from interpolateN2CH4' + write(*,*) 'Data file ',trim(dt_file),' not found.' + write(*,*) 'Check that your path to datagcm:',trim(datadir) + write(*,*) 'is correct. You can change it in callphys.def with:' + write(*,*) 'datadir = /absolute/path/to/datagcm' + write(*,*) 'Also check that the continuum data continuum_data/N2-CH4_2011.cia is there.' + call abort + else + + do iT=1,nT + + read(33,fmat1) bleh,blah,blah,nres,Ttemp + if(nS.ne.nres)then + print*,'Resolution given in file: ',trim(dt_file) + print*,'is ',nres,', which does not match nS.' + print*,'Please adjust nS value in interpolateN2CH4.F90' + stop + endif + temp_arr(iT)=Ttemp + + do k=1,nS + read(33,*) wn_arr(k),abs_arr(k,it) + end do + + end do + + endif + close(33) +!$OMP END MASTER +!$OMP BARRIER + + print*,'interpolateN2CH4: At wavenumber ',wn,' cm^-1' + print*,' temperature ',temp,' K' + print*,' N2 partial pressure ',presN2,' Pa' + print*,' and CH4 partial pressure ',presCH4,' Pa' + + endif + + call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) + +! print*,'the absorption is ',abcoef,' cm^5 molecule^-2' +! print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' + + abcoef=abcoef*losch**2*100.0*amagatN2*amagatCH4 ! convert to m^-1 + +! print*,'We have ',amagatN2,' amagats of N2' +! print*,'and ',amagatCH4,' amagats of CH4' +! print*,'So the absorption is ',abcoef,' m^-1' + + +! unlike for Rayleigh scattering, we do not currently weight by the BB function +! however our bands are normally thin, so this is no big deal. + + + return +end subroutine interpolateN2CH4 + Index: /trunk/LMDZ.TITAN/libf/phytitan/interpolateN2N2.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/interpolateN2N2.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/interpolateN2N2.F90 (revision 1648) @@ -78,5 +78,5 @@ write(*,*) 'is correct. You can change it in callphys.def with:' write(*,*) 'datadir = /absolute/path/to/datagcm' - write(*,*) 'Also check that the continuum data continuum_data/N2-N2_norm_2011.cia is there.' + write(*,*) 'Also check that the continuum data continuum_data/N2-N2_2011.cia is there.' call abort else Index: /trunk/LMDZ.TITAN/libf/phytitan/optci.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/optci.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/optci.F90 (revision 1648) @@ -1,10 +1,10 @@ subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, & QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, & - TMID,PMID,TAUGSURF,QVAR,MUVAR) + TMID,PMID,TAUGSURF,GWEIGHT) use radinc_h - use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig + use radcommon_h, only: gasi,tlimit,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig use gases_h - use comcstfi_mod, only: g, r, mugaz + use comcstfi_mod, only: g, r use callkeys_mod, only: continuum,graybody implicit none @@ -49,4 +49,17 @@ real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND) + ! Titan customisation + ! J. Vatant d'Ollone (2016) + real*8 GWEIGHT(L_NGAUSS) + real*8 DHAZE_T(L_LEVELS+1,L_NSPECTI) + real*8 DHAZES_T(L_LEVELS+1,L_NSPECTI) + real*8 SSA_T(L_LEVELS+1,L_NSPECTI) + real*8 ASF_T(L_LEVELS+1,L_NSPECTI) + real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI) + real*8 K_HAZE(L_NLAYRAD,L_NSPECTI) + + CHARACTER*2 str2 + ! ========================== + integer L, NW, NG, K, LK, IAER integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS) @@ -60,9 +73,6 @@ double precision p_cross - ! variable species mixing ratio variables - real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS) real*8 KCOEF(4) - integer NVAR(L_LEVELS) - + ! temporary variables for multiple aerosol calculation real*8 atemp @@ -73,5 +83,5 @@ !real*8 rho !! see test below - integer igas, jgas + integer igas, jgas, ilay integer interm @@ -95,10 +105,8 @@ ! if we have continuum opacities, we need dz - dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k) - U(k) = Cmk*DPR(k)*mugaz/muvar(k) ! only Cmk line in optci.F - !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present - - call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, & - LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K)) + dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K)) + U(k) = Cmk*DPR(k) ! only Cmk line in optci.F + + call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K)) do LK=1,4 @@ -120,4 +128,6 @@ do K=2,L_LEVELS + ilay = k / 2 ! int. arithmetic => gives the gcm layer index + ! continuum absorption DCONT = 0.0d0 @@ -129,9 +139,5 @@ do igas=1,ngasmx - if(gfrac(igas).eq.-1)then ! variable - p_cont = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol - else - p_cont = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k))) - endif + p_cont = dble(PMID(k)*scalep*gfrac(igas,ilay)) dtemp=0.0d0 @@ -151,9 +157,28 @@ ! then cross-interactions with other gases do jgas=1,ngasmx - p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k))) + p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) dtempc = 0.0d0 if(jgas.eq.igas_N2)then interm = indi(nw,igas,jgas) call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indi(nw,igas,jgas) = interm + endif + dtemp = dtemp + dtempc + enddo + + elseif(igas.eq.igas_CH4)then + + ! first do self-induced absorption + interm = indi(nw,igas,igas) + call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm) + indi(nw,igas,igas) = interm + + ! then cross-interactions with other gases + do jgas=1,ngasmx + p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) + dtempc = 0.0d0 + if(jgas.eq.igas_N2)then + interm = indi(nw,igas,jgas) + call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) indi(nw,igas,jgas) = interm endif @@ -185,36 +210,17 @@ DAERO=SUM(TAEROS(K,NW,1:naerkind)) +!================= Titan customisation ======================================== + call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw)) +! ============================================================================= + + do ng=1,L_NGAUSS-1 ! Now compute TAUGAS - ! Interpolate between water mixing ratios - ! WRATIO = 0.0 if the requested water amount is equal to, or outside the - ! the water data range - - if(L_REFVAR.eq.1)then ! added by RW for special no variable case - KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG) - KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG) - KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG) - KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG) - else - - KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* & - (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - & - GASI(MT(K),MP(K),NVAR(K),NW,NG)) - - KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* & - (GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - & - GASI(MT(K),MP(K)+1,NVAR(K),NW,NG)) - - KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* & - (GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - & - GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG)) - - KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* & - (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - & - GASI(MT(K)+1,MP(K),NVAR(K),NW,NG)) - - endif + KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG) + KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG) + KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG) + KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG) ! Interpolate the gaseous k-coefficients to the requested T,P values @@ -228,5 +234,6 @@ DTAUKI(K,nw,ng) = TAUGAS & + DCONT & ! For parameterized continuum absorption - + DAERO ! For aerosol absorption + + DAERO & ! For aerosol absorption + + DHAZE_T(K,NW) ! For Titan haze end do @@ -238,5 +245,6 @@ DTAUKI(K,nw,ng) = 0.d0 & + DCONT & ! For parameterized continuum absorption - + DAERO ! For aerosol absorption + + DAERO & ! For aerosol absorption + + DHAZE_T(K,NW) ! For Titan Haze end do @@ -244,7 +252,10 @@ DTAUKI(L_LEVELS+1,1:L_NSPECTI,1:L_NGAUSS)=0.d0 + + !======================================================================= ! Now the full treatment for the layers, where besides the opacity ! we need to calculate the scattering albedo and asymmetry factors + ! ====================================================================== do iaer=1,naerkind @@ -256,17 +267,34 @@ end do + ! Haze scattering DO NW=1,L_NSPECTI - DO L=1,L_NLAYRAD + DO K=2,L_LEVELS+1 + DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW) + ENDDO + ENDDO + + DO NW=1,L_NSPECTI + 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)) + 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) END DO ! L vertical loop + + !last level + L = L_NLAYRAD + K = 2*L+1 + + btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + DHAZES_T(K,NW) + END DO ! NW spectral loop +! ====================================================================== DO NW=1,L_NSPECTI NG = L_NGAUSS - DO L=1,L_NLAYRAD + DO L=1,L_NLAYRAD-1 K = 2*L+1 + DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50 @@ -278,4 +306,5 @@ GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,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) else @@ -291,18 +320,42 @@ END DO ! L vertical loop + + ! No vertical averaging on bottom layer + + L = L_NLAYRAD + K = 2*L+1 + DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + + 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) + else + WBARI(L,nw,ng) = 0.0D0 + DTAUI(L,NW,NG) = 1.0D-9 + endif + + if(btemp(L,nw) .GT. 0.0d0) then + cosbi(L,NW,NG) = atemp/btemp(L,nw) + else + cosbi(L,NW,NG) = 0.0D0 + end if ! Now the other Gauss points, if needed. DO NG=1,L_NGAUSS-1 + IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN - DO L=1,L_NLAYRAD + DO L=1,L_NLAYRAD-1 K = 2*L+1 DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50 if(DTAUI(L,NW,NG) .GT. 1.0D-9) then - - WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG) - + WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG) else WBARI(L,nw,ng) = 0.0D0 @@ -312,4 +365,18 @@ cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS) END DO ! L vertical loop + + ! No vertical averaging on bottom layer + + L = L_NLAYRAD + K = 2*L+1 + DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + if(DTAUI(L,NW,NG) .GT. 1.0D-9) then + WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG) + else + WBARI(L,nw,ng) = 0.0D0 + DTAUI(L,NW,NG) = 1.0D-9 + endif + cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS) + END IF @@ -334,16 +401,22 @@ ! ascending ray with angle theta = 0. - !open(127,file='taucum.out') - !do nw=1,L_NSPECTI - ! write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS) - !enddo - !close(127) - -! print*,'WBARI' -! print*,WBARI -! print*,'DTAUI' -! print*,DTAUI -! call abort - + +! Titan's outputs (J.V.O, 2016)=============================================== +! do l=1,L_NLAYRAD +! do nw=1,L_NSPECTI +! INT_DTAU(L,NW) = 0.0d+0 +! DO NG=1,L_NGAUSS +! INT_DTAU(L,NW)= INT_DTAU(L,NW) + dtaui(L,nw,ng)*gweight(NG) +! enddo +! enddo +! enddo + +! do nw=1,L_NSPECTI +! write(str2,'(i2.2)') nw +! call writediagfi(1,'kgi'//str2,'Gaz extinction coefficient IR band '//str2,'m-1',1,int_dtau(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1)) +! call writediagfi(1,'khi'//str2,'Haze extinction coefficient IR band '//str2,'m-1',1,k_haze(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1)) +! enddo + +! ============================================================================== return Index: /trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 (revision 1648) @@ -1,10 +1,10 @@ SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV, & QXVAER,QSVAER,GVAER,WBARV,COSBV, & - TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR) + TAURAY,TAUAERO,TMID,PMID,TAUGSURF,GWEIGHT) use radinc_h - use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig + use radcommon_h, only: gasv, tlimit, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig use gases_h - use comcstfi_mod, only: g, r, mugaz + use comcstfi_mod, only: g, r use callkeys_mod, only: continuum,graybody,callgasvis @@ -55,4 +55,17 @@ real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND) + ! Titan customisation + ! J. Vatant d'Ollone (2016) + real*8 GWEIGHT(L_NGAUSS) + real*8 DHAZE_T(L_LEVELS+1,L_NSPECTI) + real*8 DHAZES_T(L_LEVELS+1,L_NSPECTI) + real*8 SSA_T(L_LEVELS+1,L_NSPECTI) + real*8 ASF_T(L_LEVELS+1,L_NSPECTI) + real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI) + real*8 K_HAZE(L_NLAYRAD,L_NSPECTI) + + CHARACTER*2 str2 + ! ========================== + integer L, NW, NG, K, LK, IAER integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS) @@ -69,8 +82,5 @@ double precision p_cross - ! variable species mixing ratio variables - real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS) real*8 KCOEF(4) - integer NVAR(L_LEVELS) ! temporary variables for multiple aerosol calculation @@ -82,5 +92,5 @@ real*8 dz(L_LEVELS) - integer igas, jgas + integer igas, jgas, ilay integer interm @@ -107,11 +117,8 @@ ! if we have continuum opacities, we need dz - dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k) - U(k) = Cmk*DPR(k)*mugaz/muvar(k) ! only Cmk line in optci.F - !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present - - - call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, & - LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K)) + dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K)) + U(k) = Cmk*DPR(k) ! only Cmk line in optcv.F + + call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K)) do LK=1,4 @@ -133,7 +140,9 @@ end do ! levels end do - + ! we ignore K=1... do K=2,L_LEVELS + + ilay = k / 2 ! int. arithmetic => gives the gcm layer index do NW=1,L_NSPECTV @@ -153,9 +162,5 @@ do igas=1,ngasmx - if(gfrac(igas).eq.-1)then ! variable - p_cont = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol - else - p_cont = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k))) - endif + p_cont = dble(PMID(k)*scalep*gfrac(igas,ilay)) dtemp=0.0 @@ -176,5 +181,5 @@ ! then cross-interactions with other gases do jgas=1,ngasmx - p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k))) + p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) dtempc = 0.0 if(jgas.eq.igas_N2)then @@ -187,4 +192,23 @@ enddo + elseif(igas.eq.igas_CH4)then + + ! first do self-induced absorption + interm = indv(nw,igas,igas) + call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm) + indv(nw,igas,igas) = interm + + ! then cross-interactions with other gases + do jgas=1,ngasmx + p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) + dtempc = 0.0 + if(jgas.eq.igas_N2)then + interm = indv(nw,igas,jgas) + call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indv(nw,igas,jgas) = interm + endif + dtemp = dtemp + dtempc + enddo + endif @@ -197,36 +221,16 @@ endif +!================= Titan customisation ======================================== + call disr_haze(dz(k),plev(k),wnov(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw)) +! ============================================================================= + do ng=1,L_NGAUSS-1 ! Now compute TAUGAS - ! Interpolate between water mixing ratios - ! WRATIO = 0.0 if the requested water amount is equal to, or outside the - ! the water data range - - if(L_REFVAR.eq.1)then ! added by RW for special no variable case - KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG) - KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG) - KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG) - KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG) - else - - KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* & - (GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) - & - GASV(MT(K),MP(K),NVAR(K),NW,NG)) - - KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* & - (GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - & - GASV(MT(K),MP(K)+1,NVAR(K),NW,NG)) - - KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*& - (GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - & - GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG)) - - KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* & - (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - & - GASV(MT(K)+1,MP(K),NVAR(K),NW,NG)) - - endif + KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG) + KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG) + KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG) + KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG) ! Interpolate the gaseous k-coefficients to the requested T,P values @@ -240,5 +244,6 @@ DTAUKV(K,nw,ng) = TAUGAS & + TRAYAER & ! TRAYAER includes all scattering contributions - + DCONT ! For parameterized continuum aborption + + DCONT & ! For parameterized continuum aborption + + DHAZE_T(K,NW) ! For Titan haze end do @@ -248,5 +253,6 @@ NG = L_NGAUSS - DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption + DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT & ! For parameterized continuum absorption + + DHAZE_T(K,NW) ! For Titan haze do iaer=1,naerkind @@ -261,5 +267,6 @@ ! Now the full treatment for the layers, where besides the opacity ! we need to calculate the scattering albedo and asymmetry factors - + ! ====================================================================== + do iaer=1,naerkind DO NW=1,L_NSPECTV @@ -269,10 +276,24 @@ ENDDO end do + + ! Haze scattering + DO NW=1,L_NSPECTV + DO K=2,L_LEVELS+1 + DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW) + ENDDO + ENDDO + DO NW=1,L_NSPECTV 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)) - btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) + + 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) + + btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) & + + DHAZES_T(K,NW) + ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW) @@ -283,6 +304,9 @@ L = L_NLAYRAD K = 2*L+1 - atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind)) - btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + + 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) ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) @@ -292,4 +316,6 @@ END DO ! NW spectral loop +! =========================================================================================== + DO NG=1,L_NGAUSS DO NW=1,L_NSPECTV @@ -309,4 +335,5 @@ WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG) + END DO ! NW spectral loop END DO ! NG Gauss loop @@ -329,4 +356,23 @@ +! Titan's outputs (J.V.O, 2016)=============================================== +! do l=1,L_NLAYRAD +! do nw=1,L_NSPECTV +! INT_DTAU(L,NW) = 0.0d+0 +! DO NG=1,L_NGAUSS +! INT_DTAU(L,NW)= INT_DTAU(L,NW) + dtauv(L,nw,ng)*gweight(NG) +! enddo +! enddo +! enddo + +! do nw=1,L_NSPECTV +! write(str2,'(i2.2)') nw +! call writediagfi(1,'kgv'//str2,'Gaz extinction coefficient VI band '//str2,'m-1',1,int_dtau(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1)) +! call writediagfi(1,'khv'//str2,'Haze extinction coefficient VI band '//str2,'m-1',1,k_haze(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1)) +! enddo + +! ============================================================================== + + return Index: /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (revision 1648) @@ -15,5 +15,4 @@ use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind - use gases_h, only: gnom, gfrac use radcommon_h, only: sigma, glat, grav, BWNV use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe @@ -23,5 +22,5 @@ use geometry_mod, only: latitude, longitude, cell_area USE comgeomfi_h, only: totarea, totarea_planet - USE tracer_h, only: noms, mmol, radius, rho_q, qext, & + USE tracer_h, only: noms, mmol, radius, qext, & alpha_lift, alpha_devil, qextrhor, & igcm_dustbin @@ -353,6 +352,4 @@ !$OMP THREADPRIVATE(qsurf_hist) - real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW - real,allocatable,dimension(:,:,:),save :: reffrad ! Aerosol effective radius (m). By RW real,allocatable,dimension(:,:,:),save :: nueffrad ! Aerosol effective radius variance. By RW @@ -660,11 +657,11 @@ ! II.a Call correlated-k radiative transfer scheme ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - muvar(1:ngrid,1:nlayer+1)=mugaz + + call call_profilgases(nlayer) ! standard callcorrk call callcorrk(ngrid,nlayer,pq,nq,qsurf, & albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & - tsurf,fract,dist_star,aerosol,muvar, & + tsurf,fract,dist_star,aerosol, & zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, & fluxsurfabs_sw,fluxtop_lw, & Index: /trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90 (revision 1648) @@ -7,5 +7,5 @@ ! ! radcommon.h -! +!v !----------------------------------------------------------------------C ! @@ -36,6 +36,4 @@ ! TAURAY - Array (NSPECTV elements) of the pressure-independent ! part of Rayleigh scattering optical depth. -! TAURAYVAR - Array (NSPECTV elements) of the pressure-independent -! part of Rayleigh scattering optical depth for the variable gas. ! FZEROI - Fraction of zeros in the IR CO2 k-coefficients, for ! each temperature, pressure, and spectral interval @@ -63,8 +61,8 @@ REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV) !BWNV read by master in setspv - REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV), TAURAYVAR(L_NSPECTV) + REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV) !$OMP THREADPRIVATE(WNOI,DWNI,WAVEI,& !$OMP WNOV,DWNV,WAVEV,& - !$OMP STELLARF,TAURAY,TAURAYVAR) + !$OMP STELLARF,TAURAY) REAL*8 blami(L_NSPECTI+1) @@ -79,10 +77,10 @@ !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011 REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv - REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, WREFVAR, PFGASREF + REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, PFGASREF real*8 FZEROI(L_NSPECTI) real*8 FZEROV(L_NSPECTV) real*8 pgasmin, pgasmax real*8 tgasmin, tgasmax -!$OMP THREADPRIVATE(gasi,gasv,& !wrefvar,pgasref,tgasref,pfgasref read by master in sugas_corrk +!$OMP THREADPRIVATE(gasi,gasv,& !pgasref,tgasref,pfgasref read by master in sugas_corrk !$OMP FZEROI,FZEROV) !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk Index: /trunk/LMDZ.TITAN/libf/phytitan/su_gases.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/su_gases.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/su_gases.F90 (revision 1648) @@ -1,3 +1,3 @@ -subroutine su_gases +subroutine su_gases(nlayer,tracer) use gases_h @@ -5,6 +5,9 @@ implicit none - integer igas, ierr, count - + integer,intent(in) :: nlayer + logical,intent(in) :: tracer + + integer igas, ierr + !================================================================== ! @@ -17,85 +20,67 @@ ! R. Wordsworth (2011) ! Allocatable arrays by A. Spiga (2011) - ! + ! Titan's version : J. Vatant d'Ollone (2017) + ! * force igas_X labels and def gnom for N2,CH4,H2 + ! * get rid of variable species ( enrichment dimension will be set in corrk routines ) + ! !================================================================== !$OMP MASTER - ! load gas names from file 'gases.def' + + + ngasmx = 3 ! N2, CH4, H2 + + + ! load reference level and reference molar fractions from file 'gases.def' open(90,file='gases.def',status='old',form='formatted',iostat=ierr) + if (ierr.eq.0) then write(*,*) "sugases.F90: reading file gases.def" - read(90,*) - read(90,*,iostat=ierr) ngasmx + read(90,*) ! header + + ! We allocate gfrac and we set gas molar fractions for the reference level only. + ! This will be useful for the cpp_mu and rayleigh routines + ! Other gas molar fractions are now set in routine callprofilgases + + write(*,*) 'sugases.F90: allocating and reading gas molar fractions from reference level in gases.def...' + + if(.not.allocated(gfrac)) allocate(gfrac(ngasmx,nlayer)) + + read(90,*,iostat=ierr) nivref if (ierr.ne.0) then - write(*,*) "sugases.F90: error reading number of gases" + write(*,*) "sugases.F90: error reading reference level" write(*,*) " (first line of gases.def) " call abort endif - - print*,ngasmx, " gases found in gases.def. Allocating names and molar fractions..." - - if (.not.allocated(gnom)) allocate(gnom(ngasmx)) + + print*, "layer", nivref, "is reference level found in gases.def ..." + do igas=1,ngasmx - read(90,*,iostat=ierr) gnom(igas) + read(90,*,iostat=ierr) gfrac(igas,nivref) if (ierr.ne.0) then - write(*,*) 'sugases.F90: error reading gas names in gases.def...' + write(*,*) 'sugases.F90: error reading reference gas molar fractions in gases.def... aborting' call abort endif enddo !of do igas=1,ngasmx - vgas=0 - if(.not.allocated(gfrac)) allocate(gfrac(ngasmx)) - do igas=1,ngasmx - read(90,*,iostat=ierr) gfrac(igas) - if (ierr.ne.0) then - write(*,*) 'sugases.F90: error reading gas molar fractions in gases.def...' - call abort - endif - ! find variable gas (if any) - if(gfrac(igas).eq.-1.0)then - if(vgas.eq.0)then - vgas=igas - else - print*,'You seem to be choosing two variable gases' - print*,'Check that gases.def is correct' - call abort - endif - endif - - enddo !of do igas=1,ngasmx + ! We force gnom = (N2, CH4, H2) and igas_X for Titan + + igas_N2 = 1 + igas_CH4 = 2 + igas_H2 = 3 - ! assign the 'igas_X' labels - count=0 - do igas=1,ngasmx - if (trim(gnom(igas)).eq."H2_" .or. trim(gnom(igas)).eq."H2") then - igas_H2=igas - count=count+1 - elseif (trim(gnom(igas)).eq."N2_" .or. trim(gnom(igas)).eq."N2") then - igas_N2=igas - count=count+1 - elseif (trim(gnom(igas)).eq."CH4") then - igas_CH4=igas - count=count+1 - elseif (trim(gnom(igas)).eq."C2H6") then - igas_C2H6=igas - count=count+1 - elseif (trim(gnom(igas)).eq."C2H2") then - igas_C2H2=igas - count=count+1 - endif - enddo - - if(count.ne.ngasmx)then - print*,'Mismatch between ngas and number of recognised gases in sugas_corrk.F90.' - print*,'Either we haven`t managed to assign all the gases, or there are duplicates.' - print*,'Please try again.' - endif - + if (.not.allocated(gnom)) allocate(gnom(ngasmx)) + gnom(igas_N2) = "N2_" + gnom(igas_CH4) = "CH4" + gnom(igas_H2) = "H2_" + + else write(*,*) 'Cannot find required file "gases.def"' call abort endif + close(90) !$OMP END MASTER Index: /trunk/LMDZ.TITAN/libf/phytitan/sugas_corrk.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/sugas_corrk.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/sugas_corrk.F90 (revision 1648) @@ -24,5 +24,5 @@ use radcommon_h, only : tgasref,tgasmin,tgasmax use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight - use radcommon_h, only : wrefvar,WNOI,WNOV + use radcommon_h, only : WNOI,WNOV use datafile_mod, only: datadir use comcstfi_mod, only: mugaz @@ -44,5 +44,4 @@ ! ALLOCATABLE ARRAYS -- AS 12/2011 REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE,SAVE :: gasi8, gasv8 !read by master - character*20,allocatable,DIMENSION(:),SAVE :: gastype ! for check with gnom, read by master real*8 x, xi(4), yi(4), ans, p @@ -79,10 +78,4 @@ read(111,*) ngas - if(ngas.ne.ngasmx)then - print*,'Number of gases in radiative transfer data (',ngas,') does not', & - 'match that in gases.def (',ngasmx,'), exiting.' - call abort - endif - if(ngas.gt.5 .or. ngas.lt.1)then print*,ngas,' species in database [', & @@ -90,40 +83,7 @@ '], radiative code cannot handle this.' call abort - endif - - ! dynamically allocate gastype and read from Q.dat - IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( ngas ) ) - - do igas=1,ngas - read(111,*) gastype(igas) - print*,'Gas ',igas,' is ',gastype(igas) - enddo - - ! get array size, load the coefficients - open(111,file=TRIM(file_path),form='formatted') - read(111,*) L_REFVAR - IF( .NOT. ALLOCATED( wrefvar ) ) ALLOCATE( WREFVAR(L_REFVAR) ) - read(111,*) wrefvar - close(111) - - ! Check that gastype and gnom match - do igas=1,ngas - print*,'Gas ',igas,' is ',trim(gnom(igas)) - if (trim(gnom(igas)).ne.trim(gastype(igas))) then - print*,'Name of a gas in radiative transfer data (',trim(gastype(igas)),') does not ', & - 'match that in gases.def (',trim(gnom(igas)),'), exiting. You should compare ', & - 'gases.def with Q.dat in your radiative transfer directory.' - call abort - endif - enddo - print*,'Confirmed gas match in radiative transfer and gases.def!' - - ! display the values - print*,'Variable gas volume mixing ratios:' - do n=1,L_REFVAR - !print*,n,'.',wrefvar(n),' kg/kg' ! pay attention! - print*,n,'.',wrefvar(n),' mol/mol' - end do - print*,'' + endif + + L_REFVAR = 1 ! JVO 2017 : set to 1 to keep the code running until the new variable species treatment !======================================================================= @@ -628,6 +588,19 @@ enddo - endif - + elseif (igas .eq. igas_CH4) then + + ! first do self-induced absorption + dummy = -9999 + call interpolateCH4CH4(200.D+0,200.D+0,7500.D+0,testcont,.true.,dummy) + ! then cross-interactions with other gases + do jgas=1,ngasmx + if (jgas .eq. igas_N2) then + dummy = -9999 + call interpolateN2CH4(200.D+0,250.0D+0,100000.D+0,5000.D+0,testcont,.true.,dummy) + endif + enddo + + endif + enddo endif @@ -645,5 +618,4 @@ IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 ) IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref ) - IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype ) !$OMP END MASTER !$OMP BARRIER Index: /trunk/LMDZ.TITAN/libf/phytitan/tpindex.F =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/tpindex.F (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/tpindex.F (revision 1648) @@ -1,4 +1,3 @@ - subroutine tpindex(pw,tw,qvar,pref,tref,wrefvar,LCOEF,MT,MP, - & NVAR,wratio) + subroutine tpindex(pw,tw,pref,tref,LCOEF,MT,MP) !================================================================== @@ -6,5 +5,5 @@ ! Purpose ! ------- -! Interpolate K-coefficients to the given P,T and Qvar values. +! Interpolate K-coefficients to the given P,T. ! ! Notes @@ -58,8 +57,7 @@ real*8 Tref(L_NTREF) real*8 pref(L_PINT) - real*8 wrefvar(L_REFVAR) - integer MT, MP, N, M, NP, NVAR - real*8 PW, TW, Qvar, wratio + integer MT, MP, N, M, NP + real*8 PW, TW real*8 PWL, LCOEF(4), T, U @@ -143,24 +141,4 @@ LCOEF(4) = (1.0-T)*U -! Get the indicies for abundance of the varying species. There are 10 sets of -! k-coefficients with differing amounts of variable vs. constant gas. - - IF(QVAR.le.WREFVAR(1)) then - NVAR = 1 - WRATIO = 0.0D0 ! put all the weight on the first point - ELSEIF(QVAR.ge.WREFVAR(L_REFVAR)) then - NVAR = L_REFVAR-1 ! TB16 in order to not oversize NVAr when doing - !NVAR+1 - WRATIO = 1.00D0 ! put all the weight on the last point - ELSE - DO N=2,L_REFVAR - IF(QVAR.GE.WREFVAR(N-1) .and. QVAR.lt.WREFVAR(N)) then - NVAR = N-1 - WRATIO = (QVAR - WREFVAR(N-1))/(WREFVAR(N) - WREFVAR(N-1)) - GOTO 30 - END IF - END DO - END IF - 30 CONTINUE Index: /trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90 (revision 1647) +++ /trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90 (revision 1648) @@ -19,8 +19,7 @@ real varian ! Characteristic variance of log-normal distribution real r3n_q ! used to compute r0 from number and mass mixing ratio - real rho_dust ! Mars dust density (kg.m-3) real ref_r0 ! for computing reff=ref_r0*r0 (in log.n. distribution) !$OMP THREADPRIVATE(noms,mmol,radius,rho_q,qext,alpha_lift,alpha_devil,qextrhor, & - !$OMP varian,r3n_q,rho_dust,rho_ice,ref_r0) + !$OMP varian,r3n_q,ref_r0) ! tracer indexes: these are initialized in initracer and should be 0 if the