Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 772)
+++ trunk/LMDZ.GENERIC/README	(revision 773)
@@ -750,2 +750,18 @@
 - Bug fix in physiq : the size of OLR_nu is L_NSPECTI and not L_NSPECTV
 - A more robust aerosol_mod + iniaerosol : problems with ifort+parallel solved, while still OK with other compilers and seq runs.
+
+== 05/09/2012 == JL
+- Correction of the calculation of the solar longitude in tlocked case.
+    -Can now handle any prograde resonance with nres=omega_rot/omega_orb.
+    -Sun now goes westward for the standard 2:1 case, as expected.
+- In the gray case, the separation between kappa_IR and VI is now set by 
+    wave number, independently of the usual IR/VISIBLE calculation separation.
+    i.e. kappa_IR can be used in the calculation of the downward stellar flux
+      if the wavenumber in the band is low enough and vice versa.
+- In ave_stelspec, stellar flux averaging has been generalized to incorporate
+    very red/blue stellar spectra (great care must however be taken of the band
+    limit used for the corralated k distributions).
+     -Brown dwarf spectra from Allard et al. have been added.
+     -Any Black body temperature can now be used.
+
+ 
Index: trunk/LMDZ.GENERIC/libf/phystd/ave_stelspec.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/ave_stelspec.F90	(revision 772)
+++ trunk/LMDZ.GENERIC/libf/phystd/ave_stelspec.F90	(revision 773)
@@ -11,4 +11,5 @@
 !     ------- 
 !     Robin Wordsworth (2010).
+!     Generalized to very late spectral types (and Brown dwarfs) Jeremy Leconte (2012)
 !
 !     Called by
@@ -34,14 +35,13 @@
 
       integer Nfine
-      parameter(Nfine=5000)
+      integer,parameter :: Nfineband=200
       integer ifine
 
-      real lam(Nfine)
-      real stel_f(Nfine)
-      real band
+      real,allocatable :: lam(:),stel_f(:)
+      real band,lamm,lamp
       real dl
 
-      character(len=50)  :: file_id
-      character(len=100) :: file_path
+      character(len=100)  :: file_id,file_id_lam
+      character(len=200) :: file_path,file_path_lam
 
       real lam_temp
@@ -52,57 +52,84 @@
       STELLAR(:)=0.0
 
-      ! load high resolution wavenumber data
-      file_id='/stellar_spectra/lam.txt'
-      file_path=TRIM(datadir)//TRIM(file_id)
-
-      open(110,file=file_path,form='formatted',status='old',iostat=ios)
-      if (ios.ne.0) then        ! file not found
-        write(*,*) 'Error from ave_stelspec'
-        write(*,*) 'Data file ',trim(file_id),' 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 there is a ',trim(file_id),' there.'
-        call abort
-      else
-        do ifine=1,Nfine
-          read(110,*) lam(ifine)
-        enddo
-        close(110)
-      endif
-
-      dl=lam(2)-lam(1)
-
-
+      print*,'enter ave_stellspec'
       if(stelbbody)then
          tstellar=stelTbb
-         do ifine=1,Nfine
-            lam_temp=lam(ifine)
-            call blackl(dble(lam_temp*1e-6),dble(tstellar),stel_temp)
-            stel_f(ifine)=stel_temp
-         enddo
+	 Nfine=L_NSPECTV*Nfineband
+	 do band=1,L_NSPECTV
+	    lamm=10000.0/BWNV(band+1)
+	    lamp=10000.0/BWNV(band)
+	    dl=(lamp-lamm)/(Nfineband)
+	    do ifine=1,Nfineband
+	       lam_temp=lamm+(lamp-lamm)*(ifine-1.)/(Nfineband)
+               call blackl(dble(lam_temp*1e-6),dble(tstellar),stel_temp)
+	       STELLAR(band)=STELLAR(band)+stel_temp*dl
+	    enddo	    
+	 end do
+         STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
       else
          ! load high resolution stellar data
-         if(startype.eq.1)then
+         Select Case(startype)
+	   Case(1)
             file_id='/stellar_spectra/sol.txt'
             tstellar=5800.
-         elseif(startype.eq.2)then
+            file_id_lam='/stellar_spectra/lam.txt'
+            Nfine=5000
+	   Case(2)
             file_id='/stellar_spectra/gl581.txt'
             tstellar=3200.
-         elseif(startype.eq.3)then
+            file_id_lam='/stellar_spectra/lam.txt'
+            Nfine=5000
+	   Case(3)
             file_id='/stellar_spectra/adleo.txt'
             tstellar=3200.
-         elseif(startype.eq.4)then
+            file_id_lam='/stellar_spectra/lam.txt'
+            Nfine=5000
+	   Case(4)
             file_id='/stellar_spectra/gj644.txt'
             print*,'Find out tstellar before using this star!'
             call abort
-         elseif(startype.eq.5)then
+            file_id_lam='/stellar_spectra/lam.txt'
+            Nfine=5000
+	   Case(5)
             file_id='/stellar_spectra/hd128167.txt'
             tstellar=6700. ! Segura et al. (2003)
-         else
+            file_id_lam='/stellar_spectra/lam.txt'
+            Nfine=5000
+	   Case(6)
+            file_id='/stellar_spectra/BD_Teff-1600K.txt'
+            tstellar=1600. ! Segura et al. (2003)
+            file_id_lam='/stellar_spectra/lamBD.txt'
+            Nfine=5000
+	   Case(7)
+            file_id='/stellar_spectra/BD_Teff-1000K.txt'
+            tstellar=1000. ! Segura et al. (2003)
+            file_id_lam='/stellar_spectra/lamBD.txt'
+            Nfine=5000
+	   Case Default
             print*,'Error: unknown star type chosen'
             call abort
+         End Select
+
+	 allocate(lam(Nfine),stel_f(Nfine))
+
+         file_path_lam=TRIM(datadir)//TRIM(file_id_lam)
+         open(110,file=file_path_lam,form='formatted',status='old',iostat=ios)
+         if (ios.ne.0) then        ! file not found
+           write(*,*) 'Error from ave_stelspec'
+           write(*,*) 'Data file ',trim(file_id_lam),' 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 there is a ',trim(file_id_lam),' there.'
+           call abort
+         else
+           do ifine=1,Nfine
+             read(110,*) lam(ifine)
+           enddo
+           close(110)
          endif
-         
+
+
+         ! load high resolution wavenumber data
          file_path=TRIM(datadir)//TRIM(file_id)
          open(111,file=file_path,form='formatted',status='old',iostat=ios)
@@ -121,24 +148,27 @@
            close(111)
          endif
+	 
+         ! sum data by band
+         band=1
+	 Do while(lam(1).lt. real(10000.0/BWNV(band+1)))
+	    if (band.gt.L_NSPECTV) exit
+            band=band+1
+	 enddo
+	 dl=lam(2)-lam(1)
+         STELLAR(band)=STELLAR(band)+stel_f(1)*dl
+         do ifine = 2,Nfine
+            if(lam(ifine) .gt. real(10000.0/BWNV(band)))then
+               band=band-1
+            endif
+            if(band .lt. 1) exit
+	    dl=lam(ifine)-lam(ifine-1)
+            STELLAR(band)=STELLAR(band)+stel_f(ifine)*dl
+         end do
+	       
+	 
+         STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
+         deallocate(lam,stel_f)
+         
       endif
 
-
-      ! sum data by band
-      band=1
-      do ifine = 1,Nfine
-
-         if(lam(Nfine-ifine+1) .lt. real(10000.0/BWNV(band+1)))then
-            band=band+1
-         endif
-         if(band .gt. L_NSPECTV)then
-            goto 9999 ! ok, ok, I know they're evil
-         endif
-         STELLAR(band)=STELLAR(band)+stel_f(Nfine-ifine+1)*dl
-
-      end do
-
-9999  continue
-      STELLAR=STELLAR/sum(STELLAR)
-
-
       end subroutine ave_stelspec
Index: trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 772)
+++ trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 773)
@@ -194,6 +194,5 @@
                hice(ig)        = max(hice(ig),0.0)
                hice(ig)        = min(hice(ig),maxicethick)
-               pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT* &
-                   rhowater/pcapcal(ig)/ptimestep              
+               pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep              
                albedo(ig)      = albedoice
 
Index: trunk/LMDZ.GENERIC/libf/phystd/largescale.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 772)
+++ trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 773)
@@ -91,5 +91,8 @@
       DO i = 1, ngridmx
 
-         if(zt(i).le.15.) zt(i)=15.   ! check too low temperatures
+         if(zt(i).le.15.) then
+	    print*,'in lsc',i,zt(i)
+	    zt(i)=15.   ! check too low temperatures
+         endif
 !         call watersat(zt(i),pplay(i,k),zqs(i))
          call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
Index: trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 772)
+++ trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 773)
@@ -50,5 +50,5 @@
 !     Local variables
       INTEGER i, k, iq, nn
-      INTEGER, PARAMETER :: niter = 6
+      INTEGER, PARAMETER :: niter = 1
       INTEGER k1, k1p, k2, k2p
       LOGICAL itest(ngridmx)
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 772)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 773)
@@ -642,6 +642,9 @@
            if (tlocked) then
               ztim1=SIN(declin)
-              ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres)
-              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
+!              ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres)
+!              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
+! JL12 corrects tidally resonant cases. nres=omega_rot/omega_orb
+              ztim2=COS(declin)*COS(2.*pi*(zday/year_day)*(nres-1.))
+              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day)*(nres-1.))
 
               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
@@ -1017,9 +1020,5 @@
 !             Moist convection
 !             ----------------
-               ! Re-evaporate cloud water/ice
-!	       dqevap1(1:ngridmx,1:nlayermx)=0.
-!	       dtevap1(1:ngridmx,1:nlayermx)=0.
-!               call evap(ptimestep,pt,pq,pdq,pdt,dqevap1,dtevap1,qevap,tevap)
-               
+
 	       dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0.
 	       dtmoist(1:ngridmx,1:nlayermx)=0.
@@ -1032,14 +1031,4 @@
                           +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice)
                pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx)
-
-!               do l=1,nlayermx
-!	         do ig=1,ngridmx
-!	          if(isnan(dqmoist(ig,l,igcm_h2o_ice))) print*,'Nan in dqmoist ice, abort,ig,l',ig,l
-!	          if(isnan(dqmoist(ig,l,igcm_h2o_vap))) print*,'Nan in dqmoist vap, abort,ig,l',ig,l
-!	          if(isnan(dtmoist(ig,l))) print*,'Nan in dtmoist, abort,ig,l',ig,l		  
-!                  if(isnan(dqmoist(ig,l,igcm_h2o_ice)).or.isnan(dqmoist(ig,l,igcm_h2o_vap)).or.isnan(dtmoist(ig,l))) STOP
-!	         enddo
-!	       enddo
-!	       print*,'after moist',dqmoist(185,1,igcm_h2o_ice),dqmoist(185,1,igcm_h2o_vap),dtmoist(185,1)*86400.
 
                !-------------------------
@@ -1051,5 +1040,6 @@
                   enddo
                   print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
-                  print*,'In moistadj MAX atmospheric energy change   =',MAXVAL(pdt(:,:))*86400.,'K/day'
+                  print*,'In moistadj MAX atmospheric energy change   =',MAXVAL(dtmoist(:,:))*ptimestep,'K/step'
+                  print*,'In moistadj MIN atmospheric energy change   =',MINVAL(dtmoist(:,:))*ptimestep,'K/step'
 		  
 		! test energy conservation
Index: trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 772)
+++ trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 773)
@@ -299,4 +299,5 @@
 !!! - physical constants: nevermind, things are done allright below
 !!! - physical frequency: nevermind, in inifis this is a simple print
+      cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
       CALL inifis(1,llm,day0,daysec,dtphys,
      .            lati,long,area,rad,g,r,cpp)
Index: trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 772)
+++ trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 773)
@@ -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
+      use radcommon_h, only : wrefvar,WNOI,WNOV
       use datafile_mod, only: datadir
 
@@ -49,5 +49,6 @@
 
       real*8 x, xi(4), yi(4), ans, p
-      real kappa_IR, kappa_VI
+!     For gray case (JL12)
+      real kappa_IR, kappa_VI, IR_VI_wnlimit, nVI_limit,nIR_limit
 
       integer ngas, igas, jgas
@@ -274,6 +275,44 @@
 !     Get gaseous k-coefficients and interpolate onto finer pressure grid
 
+      if (graybody) then
+!        constant absorption coefficient in visible
+         write(*,*)"graybody: constant absorption coefficient in visible:"
+         kappa_VI=-100000.
+         call getin("kappa_VI",kappa_VI)
+         write(*,*)" kappa_VI = ",kappa_VI
+	 kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule         
+      
+!        constant absorption coefficient in IR
+         write(*,*)"graybody: constant absorption coefficient in InfraRed:"
+         kappa_IR=-100000.
+         call getin("kappa_IR",kappa_IR)
+         write(*,*)" kappa_IR = ",kappa_IR	 
+         kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule 
+	       
+!        wavelength used to separate IR from VI
+         IR_VI_wnlimit=3000.
+         write(*,*)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
+	 
+	 nVI_limit=L_NSPECTV
+	 Do nw=1,L_NSPECTV
+	    if (WNOV(nw).gt.IR_VI_wnlimit) then
+	       nVI_limit=nw-1
+	       exit
+	    endif
+	 End do
+	 nIR_limit=L_NSPECTI
+	 Do nw=1,L_NSPECTI
+	    if (WNOI(nw).gt.IR_VI_wnlimit) then
+	       nIR_limit=nw-1
+	       exit
+	    endif
+	 End do
+         write(*,*)"graybody: Visible / Infrared separation set at band ",nIR_limit,nVI_limit
+	      
+      End if
+
+
       ! VISIBLE
-      if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not.graybody)) then
+      if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then
          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat' 
          file_path=TRIM(datadir)//TRIM(file_id)
@@ -293,20 +332,19 @@
 
       else if (callgasvis.and.graybody) then
-         ! constant absorption coefficient in visible
-         write(*,*)"constant absorption coefficient in visible:"
-         kappa_VI = -100000.
-         call getin("kappa_VI",kappa_VI)
-         write(*,*)" kappa_VI = ",kappa_VI
-	 kappa_VI = kappa_VI * 1.e4 * mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule         
-         gasv8(:,:,:,:,:) = kappa_VI
-
+         if(nVI_limit.eq.0) then
+	    gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_VI
+         else if (nVI_limit.eq.L_NSPECTV) then
+	    gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_IR
+	 else
+	    gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)=kappa_IR
+	    gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)=kappa_VI
+	 end if
       else
          print*,'Visible corrk gaseous absorption is set to zero.'
-         gasv8(:,:,:,:,:) = 0.0
-
+         gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
       endif
 
       ! INFRA-RED
-      if ((.not.TRIM(corrkdir).eq.'null').and.(.not.graybody)) then
+      if ((.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then
          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat' 
          file_path=TRIM(datadir)//TRIM(file_id)
@@ -364,16 +402,15 @@
 
       else if (graybody) then
-         ! constant absorption coefficient in IR
-         write(*,*)"constant absorption coefficient in InfraRed:"
-         kappa_IR = -100000.
-         call getin("kappa_IR",kappa_IR)
-         write(*,*)" kappa_IR = ",kappa_IR
-         kappa_IR = kappa_IR * 1.e4 * mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule       
-         gasi8(:,:,:,:,:) = kappa_IR
-
+         if(nIR_limit.eq.0) then
+	    gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_VI
+         else if (nIR_limit.eq.L_NSPECTI) then
+	    gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_IR
+	 else
+	    gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)=kappa_IR
+	    gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)=kappa_VI
+         end if
       else
          print*,'Infrared corrk gaseous absorption is set to zero.'
-         gasi8(:,:,:,:,:) = 0.0
-
+         gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
       endif
 
@@ -389,6 +426,5 @@
                   do nw=1,L_NSPECTV
                      if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
-                        gasv8(nt,np,nh,nw,ng) = &
-                            log10(gasv8(nt,np,nh,nw,ng))
+                        gasv8(nt,np,nh,nw,ng) = log10(gasv8(nt,np,nh,nw,ng))
                      else
                         gasv8(nt,np,nh,nw,ng) = -200.0
@@ -398,6 +434,5 @@
                   do nw=1,L_NSPECTI
                      if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
-                        gasi8(nt,np,nh,nw,ng) = &
-                            log10(gasi8(nt,np,nh,nw,ng))
+                        gasi8(nt,np,nh,nw,ng) = log10(gasi8(nt,np,nh,nw,ng))
                      else
                         gasi8(nt,np,nh,nw,ng) = -200.0
