Index: trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90 (revision 873) +++ trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90 (revision 873) @@ -0,0 +1,84 @@ + subroutine bilinearbig(x_arr,y_arr,f2d_arr,x_in,y_in,f,a) + +! Necessary for interpolation of continuum data +! optimized by A. Spiga 01/2013 + + implicit none + + integer nX,nY,i,j,a,b + parameter(nX=2428) + parameter(nY=10) + + real*8 x_in,y_in,x1,x2,y1,y2 + real*8 f,f11,f12,f21,f22,fA,fB + real*8 x_arr(nX) + real*8 y_arr(nY) + real*8 f2d_arr(nX,nY) + real*8,save :: x,y + + integer strlen + character*100 label + label='subroutine bilinear' + + x=x_in + y=y_in + + + !! AS: important to optimize here because the array is quite large + !! ... and actually calculations only need to be done once + !! --> Case 1 : we have not calculated yet + if ( a == -9999) then + !1st check we're within the wavenumber range + if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then + f=0.0D+0 + a=-1 + else + i=1 + x2=x_arr(i) + do while ( x2 .le. x ) + x1=x2 + i=i+1 + x2=x_arr(i) + a=i-1 + end do + endif + !! --> case 2 : we already saw we are out of wavenumber range + else if ( a == -1) then + f=0.0D+0 + return + !! --> case 3 : we already determined a -- so we just proceed + else + x1=x_arr(a) + x2=x_arr(a+1) + endif + +! ... and for y within the temperature range + if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then + write(*,*) 'Warning from bilinearH2H2:' + write(*,*) 'Outside continuum temperature range!' + if(y.lt.y_arr(1))then + y=y_arr(1)+0.01 + endif + if(y.gt.y_arr(nY))then + y=y_arr(nY)-0.01 + endif + else + j=1 + y2=y_arr(j) + do while ( y2 .le. y ) + y1=y2 + j=j+1 + y2=y_arr(j) + b=j-1 + end do + endif + + f11=f2d_arr(a,b) + f21=f2d_arr(a+1,b) + f12=f2d_arr(a,b+1) + f22=f2d_arr(a+1,b+1) + + call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) + + return + end subroutine bilinearbig Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys.h =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/callkeys.h (revision 869) +++ trunk/LMDZ.GENERIC/libf/phystd/callkeys.h (revision 873) @@ -8,5 +8,5 @@ & , iaervar,iddist,topdustref,callstats,calleofdump & & , enertest & - & , callgasvis,Continuum,H2Ocont_simple,graybody & + & , callgasvis,continuum,H2Ocont_simple,graybody & & , Nmix_co2, radfixed, dusttau & & , meanOLR, specOLR & @@ -35,5 +35,5 @@ & , season,diurnal,tlocked,lwrite & & , callstats,calleofdump & - & , callgasvis,Continuum,H2Ocont_simple,graybody + & , callgasvis,continuum,H2Ocont_simple,graybody logical enertest Index: trunk/LMDZ.GENERIC/libf/phystd/inifis.F =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/inifis.F (revision 869) +++ trunk/LMDZ.GENERIC/libf/phystd/inifis.F (revision 873) @@ -203,7 +203,7 @@ write(*,*) "call continuum opacities in radiative transfer ?", & "(matters only if callrad=T)" - Continuum=.true. ! default value - call getin("Continuum",Continuum) - write(*,*) " Continuum = ",Continuum + continuum=.true. ! default value + call getin("continuum",continuum) + write(*,*) " continuum = ",continuum write(*,*) "use analytic function for H2O continuum ?" Index: trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90 (revision 869) +++ trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90 (revision 873) @@ -1,3 +1,3 @@ - subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall) + subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,a) !================================================================== @@ -15,4 +15,5 @@ use datafile_mod, only: datadir + implicit none @@ -52,4 +53,6 @@ double precision blah, Ttemp integer nres + + integer a if(temp.gt.400)then @@ -104,124 +107,19 @@ print*,' pressure ',pres,' Pa' - call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef) + endif - print*,'the absorption is ',abcoef,' cm^5 molecule^-2' - print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' + call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a) + + !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 - print*,'We have ',amagat,' amagats of H2' - print*,'So the absorption is ',abcoef,' m^-1' - - !open(117,file='T_array.dat') - !do iT=1,nT - ! write(117,*), temp_arr(iT) - !end do - !close(117) - - !open(115,file='nu_array.dat') - !do k=1,nS - ! write(115,*), wn_arr(k) - !end do - !close(115) - - !open(113,file='abs_array.dat') - !do iT=1,nT - ! do k=1,nS - ! write(113,*), abs_arr(k,iT)*losch**2 - ! end do - !end do - !close(113) - - else - - call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef) - abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1 + !print*,'We have ',amagat,' amagats of H2' + !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. - endif - return end subroutine interpolateH2H2 - - -!------------------------------------------------------------------------- - subroutine bilinearH2H2(x_arr,y_arr,f2d_arr,x_in,y_in,f) -! Necessary for interpolation of continuum data - - implicit none - - integer nX,nY,i,j,a,b - parameter(nX=2428) - parameter(nY=10) - - real*8 x_in,y_in,x,y,x1,x2,y1,y2 - real*8 f,f11,f12,f21,f22,fA,fB - real*8 x_arr(nX) - real*8 y_arr(nY) - real*8 f2d_arr(nX,nY) - - integer strlen - character*100 label - label='subroutine bilinear' - - x=x_in - y=y_in - - -! 1st check we're within the wavenumber range - if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then - f=0.0D+0 - return - else - -! in the x (wavenumber) direction 1st - i=1 - 10 if (i.lt.(nX+1)) then - if (x_arr(i).gt.x) then - x1=x_arr(i-1) - x2=x_arr(i) - a=i-1 - i=9999 - endif - i=i+1 - goto 10 - endif - endif - - if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then - write(*,*) 'Warning from bilinearH2H2:' - write(*,*) 'Outside continuum temperature range!' - if(y.lt.y_arr(1))then - y=y_arr(1)+0.01 - endif - if(y.gt.y_arr(nY))then - y=y_arr(nY)-0.01 - endif - else - -! in the y (temperature) direction 2nd - j=1 - 20 if (j.lt.(nY+1)) then - if (y_arr(j).gt.y) then - y1=y_arr(j-1) - y2=y_arr(j) - b=j-1 - j=9999 - endif - j=j+1 - goto 20 - endif - endif - - f11=f2d_arr(a,b) - f21=f2d_arr(a+1,b) - f12=f2d_arr(a,b+1) - f22=f2d_arr(a+1,b+1) - - call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) - - return - end subroutine bilinearH2H2 Index: trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90 (revision 869) +++ trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90 (revision 873) @@ -1,3 +1,3 @@ - subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall) + subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall,a) !================================================================== @@ -15,4 +15,5 @@ use datafile_mod, only: datadir + implicit none @@ -55,5 +56,6 @@ integer nres - + integer a + if(temp.gt.400)then print*,'Your temperatures are too high for this H2-He CIA dataset.' @@ -108,105 +110,20 @@ print*,' and He partial pressure ',presHe,' Pa' - call bilinearH2He(wn_arr,temp_arr,abs_arr,wn,temp,abcoef) + endif - print*,'the absorption is ',abcoef,' cm^5 molecule^-2' - print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' + call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a) + + !print*,'the absorption is ',abcoef,' cm^5 molecule^-2' + !print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' abcoef=abcoef*losch**2*100.0*amagatH2*amagatHe ! convert to m^-1 - print*,'We have ',amagatH2,' amagats of H2' - print*,'and ',amagatHe,' amagats of He' - print*,'So the absorption is ',abcoef,' m^-1' - - else - - call bilinearH2He(wn_arr,temp_arr,abs_arr,wn,temp,abcoef) - abcoef=abcoef*losch**2*100.0*amagatH2*amagatHe ! convert to m^-1 + !print*,'We have ',amagatH2,' amagats of H2' + !print*,'and ',amagatHe,' amagats of He' + !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. - endif - return end subroutine interpolateH2He - - -!------------------------------------------------------------------------- - subroutine bilinearH2He(x_arr,y_arr,f2d_arr,x_in,y_in,f) -! Necessary for interpolation of continuum data - - implicit none - - integer nX,nY,i,j,a,b - parameter(nX=2428) - parameter(nY=10) - - real*8 x_in,y_in,x,y,x1,x2,y1,y2 - real*8 f,f11,f12,f21,f22,fA,fB - real*8 x_arr(nX) - real*8 y_arr(nY) - real*8 f2d_arr(nX,nY) - - integer strlen - character*100 label - label='subroutine bilinear' - - x=x_in - y=y_in - -! 1st check we're within the wavenumber range - if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then - f=0.0D+0 - return - else - -! in the x (wavenumber) direction 1st - i=1 - 10 if (i.lt.(nX+1)) then - if (x_arr(i).gt.x) then - x1=x_arr(i-1) - x2=x_arr(i) - a=i-1 - i=9999 - endif - i=i+1 - goto 10 - endif - endif - - if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then - write(*,*) 'Warning from bilinearH2He:' - write(*,*) 'Outside continuum temperature range!' - if(y.lt.y_arr(1))then - y=y_arr(1)+0.01 - endif - if(y.gt.y_arr(nY))then - y=y_arr(nY)-0.01 - endif - else - -! in the y (temperature) direction 2nd - j=1 - 20 if (j.lt.(nY+1)) then - if (y_arr(j).gt.y) then - y1=y_arr(j-1) - y2=y_arr(j) - b=j-1 - j=9999 - endif - j=j+1 - goto 20 - endif - endif - - f11=f2d_arr(a,b) - f21=f2d_arr(a+1,b) - f12=f2d_arr(a,b+1) - f22=f2d_arr(a+1,b+1) - - call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) - - - return - end subroutine bilinearH2He Index: trunk/LMDZ.GENERIC/libf/phystd/optci.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/optci.F90 (revision 869) +++ trunk/LMDZ.GENERIC/libf/phystd/optci.F90 (revision 873) @@ -4,5 +4,5 @@ use radinc_h - use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep + use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi use gases_h implicit none @@ -70,7 +70,10 @@ ! variables for k in units m^-1 - real*8 rho, dz(L_LEVELS) + real*8 dz(L_LEVELS) + !real*8 rho !! see test below integer igas, jgas + + integer interm !--- Kasting's CIA ---------------------------------------- @@ -87,4 +90,10 @@ !---------------------------------------------------------- + !! AS: to save time in computing continuum (see bilinearbig) + IF (.not.ALLOCATED(indi)) THEN + ALLOCATE(indi(L_NSPECTI,ngasmx)) + indi = -9999 ! this initial value means "to be calculated" + ENDIF + !======================================================================= ! Determine the total gas opacity throughout the column, for each @@ -133,9 +142,10 @@ do K=2,L_LEVELS - do nw=1,L_NSPECTI + + do NW=1,L_NSPECTI DCONT = 0.0 ! continuum absorption - if(Continuum.and.(.not.graybody))then + if(continuum.and.(.not.graybody))then ! include continua if necessary wn_cont = dble(wnoi(nw)) @@ -157,5 +167,7 @@ ! first do self-induced absorption - call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.) + interm = indi(nw,igas) + call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) + indi(nw,igas) = interm ! then cross-interactions with other gases @@ -166,5 +178,7 @@ call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.) elseif(jgas.eq.igas_He)then - call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.) + interm = indi(nw,jgas) + call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indi(nw,jgas) = interm endif dtemp = dtemp + dtempc @@ -228,5 +242,5 @@ 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)+1,NW,NG) - & GASI(MT(K),MP(K),NVAR(K),NW,NG)) @@ -242,4 +256,5 @@ (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - & GASI(MT(K)+1,MP(K),NVAR(K),NW,NG)) + endif @@ -252,5 +267,6 @@ TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT - DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption + DTAUKI(K,nw,ng) = TAUGAS & + + DCONT ! For parameterized continuum absorption do iaer=1,naerkind @@ -278,11 +294,11 @@ ! we need to calculate the scattering albedo and asymmetry factors + do iaer=1,naerkind DO NW=1,L_NSPECTI DO K=2,L_LEVELS+1 - do iaer=1,naerkind TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) - end do ENDDO ENDDO + end do DO NW=1,L_NSPECTI Index: trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 (revision 869) +++ trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 (revision 873) @@ -4,5 +4,5 @@ use radinc_h - use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep + use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv use gases_h @@ -35,6 +35,7 @@ !------------------------------------------------------------------- +#include "comcstfi.h" #include "callkeys.h" -#include "comcstfi.h" + real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) @@ -46,26 +47,30 @@ real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) - real*8 TAURAY(L_NSPECTV) ! for aerosols - real*8 QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND) - real*8 QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND) - real*8 GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND) - real*8 TAUAERO(L_LEVELS+1,NAERKIND) - real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND) - real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND) - - integer L, NW, NG, K, NG1(L_NSPECTV), LK, IAER + real*8 QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND) + real*8 QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND) + real*8 GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND) + real*8 TAUAERO(L_LEVELS+1,NAERKIND) + real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND) + real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND) + + integer L, NW, NG, K, LK, IAER integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS) real*8 ANS, TAUGAS + real*8 TAURAY(L_NSPECTV) real*8 TRAY(L_LEVELS,L_NSPECTV) + real*8 TRAYAER real*8 DPR(L_LEVELS), U(L_LEVELS) real*8 LCOEF(4), LKCOEF(L_LEVELS,4) - real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER + real*8 taugsurf(L_NSPECTV,L_NGAUSS-1) + real*8 DCONT + double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc + double precision p_cross ! variable species mixing ratio variables - real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS) - real*8 KCOEF(4) + real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS) + real*8 KCOEF(4) integer NVAR(L_LEVELS) @@ -74,9 +79,15 @@ ! variables for k in units m^-1 - double precision wn_cont, p_cont, p_air, T_cont, dtemp - double precision p_cross - real*8 dz(L_LEVELS), DCONT + real*8 dz(L_LEVELS) integer igas, jgas + + integer interm + + !! AS: to save time in computing continuum (see bilinearbig) + IF (.not.ALLOCATED(indv)) THEN + ALLOCATE(indv(L_NSPECTV,ngasmx)) + indv = -9999 ! this initial value means "to be calculated" + ENDIF !======================================================================= @@ -95,5 +106,5 @@ ! if we have continuum opacities, we need dz if(kastprof)then - dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K)) + dz(k) = dpr(k)*(1000.0*8.3145/muvar(k))*TMID(K)/(g*PMID(K)) U(k) = (Cmk*mugaz/(muvar(k)))*DPR(k) else @@ -109,21 +120,20 @@ end do + DO NW=1,L_NSPECTV + do iaer=1,naerkind + TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER) + end do TRAY(K,NW) = TAURAY(NW) * DPR(K) - - do iaer=1,naerkind - TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER) - end do - END DO - end do - - ! TRAYAER is Tau RAYleigh scattering, plus AERosol opacity + end do ! levels ! we ignore K=1... do K=2,L_LEVELS + do NW=1,L_NSPECTV TRAYAER = TRAY(K,NW) + ! TRAYAER is Tau RAYleigh scattering, plus AERosol opacity do iaer=1,naerkind TRAYAER = TRAYAER + TAEROS(K,NW,IAER) @@ -132,5 +142,5 @@ DCONT = 0.0 ! continuum absorption - if(callgasvis.and.continuum.and.(.not.graybody))then + if(continuum.and.(.not.graybody).and.callgasvis)then ! include continua if necessary wn_cont = dble(wnov(nw)) @@ -153,15 +163,21 @@ ! first do self-induced absorption - call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.) + interm = indv(nw,igas) + call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) + indv(nw,igas) = interm ! then cross-interactions with other gases do jgas=1,ngasmx p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k))) - if(jgas.eq.igas_N2)then - call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.) + dtempc = 0.0 + if(jgas.eq.igas_N2)then + call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.) ! should be irrelevant in the visible elseif(jgas.eq.igas_He)then - call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.) + interm = indv(nw,jgas) + call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indv(nw,jgas) = interm endif + dtemp = dtemp + dtempc enddo @@ -181,16 +197,17 @@ enddo - DCONT = DCONT*dz(k) + DCONT = DCONT*dz(k) + endif - 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 + 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) @@ -198,4 +215,5 @@ 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) - & @@ -213,29 +231,29 @@ (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - & GASV(MT(K)+1,MP(K),NVAR(K),NW,NG)) + endif - ! Interpolate the gaseous k-coefficients to the requested T,P values - - ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) + & + ! Interpolate the gaseous k-coefficients to the requested T,P values + + ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) + & LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4) - TAUGAS = U(k)*ANS - - !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + TAUGAS = U(k)*ANS + TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT - DTAUKV(K,nw,ng) = TAUGAS + TRAYAER & ! TRAYAER includes all scattering contributions - + DCONT ! for continuum absorption + DTAUKV(K,nw,ng) = TAUGAS & + + TRAYAER & ! TRAYAER includes all scattering contributions + + DCONT ! For parameterized continuum aborption end do - - ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS), - ! which holds continuum opacity only - - NG = L_NGAUSS + ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS), + ! which holds continuum opacity only + + NG = L_NGAUSS DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption + do iaer=1,naerkind - DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER) - ! & + DCONT ! For parameterized continuum absorption + DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER) end do ! a bug was here! @@ -248,37 +266,33 @@ ! we need to calculate the scattering albedo and asymmetry factors + do iaer=1,naerkind DO NW=1,L_NSPECTV - DO K=2,L_LEVELS - do iaer=1,naerkind + DO K=2,L_LEVELS ! AS: shouldn't this be L_LEVELS+1 ? (see optci) TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) - end do ENDDO ENDDO - + end do DO NW=1,L_NSPECTV - DO NG=1,L_NGAUSS - DO L=1,L_NLAYRAD-1 - K = 2*L+1 - - DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)+DTAUKV(K+1,NW,NG) - - atemp=0. - btemp=TRAY(K,NW) + TRAY(K+1,NW) - ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) + DO NG=1,L_NGAUSS + DO L=1,L_NLAYRAD-1 + + K = 2*L+1 + DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG) + + atemp = 0. + btemp = TRAY(K,NW) + TRAY(K+1,NW) + ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) do iaer=1,naerkind atemp = atemp + & GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + & GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER) - btemp = btemp + & - TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER) - ctemp = ctemp + & - TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER) + btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER) + ctemp = ctemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER) end do - + WBARV(L,nw,ng) = ctemp / DTAUV(L,nw,ng) COSBV(L,NW,NG) = atemp/btemp - WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng) - - END DO + + END DO ! L vertical loop ! No vertical averaging on bottom layer @@ -299,12 +313,10 @@ WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng) - END DO ! NG gauss point loop + END DO ! NG Gauss loop END DO ! NW spectral loop - - ! Total extinction optical depths - DO NW=1,L_NSPECTV + DO NW=1,L_NSPECTV DO NG=1,L_NGAUSS ! full gauss loop TAUV(1,NW,NG)=0.0D0 @@ -321,4 +333,6 @@ - RETURN -END SUBROUTINE OPTCV + return + + +end subroutine optcv Index: trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90 (revision 869) +++ trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90 (revision 873) @@ -66,4 +66,8 @@ REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90 + !! AS: introduced to avoid doing same computations again for continuum + INTEGER, DIMENSION(:,:), ALLOCATABLE :: indi + INTEGER, DIMENSION(:,:), ALLOCATABLE :: indv + !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011 REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv Index: trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90 (revision 869) +++ trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90 (revision 873) @@ -55,4 +55,6 @@ double precision testcont ! for continuum absorption initialisation + + integer :: dummy !======================================================================= @@ -596,5 +598,6 @@ ! first do self-induced absorption - call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.) + dummy = -9999 + call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.,dummy) ! then cross-interactions with other gases do jgas=1,ngasmx @@ -602,5 +605,6 @@ call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.) elseif (jgas .eq. igas_He) then - call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.) + dummy = -9999 + call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy) endif enddo