Changeset 773 for trunk/LMDZ.GENERIC


Ignore:
Timestamp:
Sep 5, 2012, 8:32:32 PM (12 years ago)
Author:
jleconte
Message:

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.


Location:
trunk/LMDZ.GENERIC
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r751 r773  
    750750- Bug fix in physiq : the size of OLR_nu is L_NSPECTI and not L_NSPECTV
    751751- A more robust aerosol_mod + iniaerosol : problems with ifort+parallel solved, while still OK with other compilers and seq runs.
     752
     753== 05/09/2012 == JL
     754- Correction of the calculation of the solar longitude in tlocked case.
     755    -Can now handle any prograde resonance with nres=omega_rot/omega_orb.
     756    -Sun now goes westward for the standard 2:1 case, as expected.
     757- In the gray case, the separation between kappa_IR and VI is now set by
     758    wave number, independently of the usual IR/VISIBLE calculation separation.
     759    i.e. kappa_IR can be used in the calculation of the downward stellar flux
     760      if the wavenumber in the band is low enough and vice versa.
     761- In ave_stelspec, stellar flux averaging has been generalized to incorporate
     762    very red/blue stellar spectra (great care must however be taken of the band
     763    limit used for the corralated k distributions).
     764     -Brown dwarf spectra from Allard et al. have been added.
     765     -Any Black body temperature can now be used.
     766
     767 

  • trunk/LMDZ.GENERIC/libf/phystd/ave_stelspec.F90

    r374 r773  
    1111!     -------
    1212!     Robin Wordsworth (2010).
     13!     Generalized to very late spectral types (and Brown dwarfs) Jeremy Leconte (2012)
    1314!
    1415!     Called by
    … …  
    3435
    3536      integer Nfine
    36       parameter(Nfine=5000)
     37      integer,parameter :: Nfineband=200
    3738      integer ifine
    3839
    39       real lam(Nfine)
    40       real stel_f(Nfine)
    41       real band
     40      real,allocatable :: lam(:),stel_f(:)
     41      real band,lamm,lamp
    4242      real dl
    4343
    44       character(len=50)  :: file_id
    45       character(len=100) :: file_path
     44      character(len=100)  :: file_id,file_id_lam
     45      character(len=200) :: file_path,file_path_lam
    4646
    4747      real lam_temp
    … …  
    5252      STELLAR(:)=0.0
    5353
    54       ! load high resolution wavenumber data
    55       file_id='/stellar_spectra/lam.txt'
    56       file_path=TRIM(datadir)//TRIM(file_id)
    57 
    58       open(110,file=file_path,form='formatted',status='old',iostat=ios)
    59       if (ios.ne.0) then        ! file not found
    60         write(*,*) 'Error from ave_stelspec'
    61         write(*,*) 'Data file ',trim(file_id),' not found.'
    62         write(*,*)'Check that your path to datagcm:',trim(datadir)
    63         write(*,*)' is correct. You can change it in callphys.def with:'
    64         write(*,*)' datadir = /absolute/path/to/datagcm'
    65         write(*,*)' Also check that there is a ',trim(file_id),' there.'
    66         call abort
    67       else
    68         do ifine=1,Nfine
    69           read(110,*) lam(ifine)
    70         enddo
    71         close(110)
    72       endif
    73 
    74       dl=lam(2)-lam(1)
    75 
    76 
     54      print*,'enter ave_stellspec'
    7755      if(stelbbody)then
    7856         tstellar=stelTbb
    79          do ifine=1,Nfine
    80             lam_temp=lam(ifine)
    81             call blackl(dble(lam_temp*1e-6),dble(tstellar),stel_temp)
    82             stel_f(ifine)=stel_temp
    83          enddo
     57         Nfine=L_NSPECTV*Nfineband
     58         do band=1,L_NSPECTV
     59            lamm=10000.0/BWNV(band+1)
     60            lamp=10000.0/BWNV(band)
     61            dl=(lamp-lamm)/(Nfineband)
     62            do ifine=1,Nfineband
     63               lam_temp=lamm+(lamp-lamm)*(ifine-1.)/(Nfineband)
     64               call blackl(dble(lam_temp*1e-6),dble(tstellar),stel_temp)
     65               STELLAR(band)=STELLAR(band)+stel_temp*dl
     66            enddo           
     67         end do
     68         STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
    8469      else
    8570         ! load high resolution stellar data
    86          if(startype.eq.1)then
     71         Select Case(startype)
     72           Case(1)
    8773            file_id='/stellar_spectra/sol.txt'
    8874            tstellar=5800.
    89          elseif(startype.eq.2)then
     75            file_id_lam='/stellar_spectra/lam.txt'
     76            Nfine=5000
     77           Case(2)
    9078            file_id='/stellar_spectra/gl581.txt'
    9179            tstellar=3200.
    92          elseif(startype.eq.3)then
     80            file_id_lam='/stellar_spectra/lam.txt'
     81            Nfine=5000
     82           Case(3)
    9383            file_id='/stellar_spectra/adleo.txt'
    9484            tstellar=3200.
    95          elseif(startype.eq.4)then
     85            file_id_lam='/stellar_spectra/lam.txt'
     86            Nfine=5000
     87           Case(4)
    9688            file_id='/stellar_spectra/gj644.txt'
    9789            print*,'Find out tstellar before using this star!'
    9890            call abort
    99          elseif(startype.eq.5)then
     91            file_id_lam='/stellar_spectra/lam.txt'
     92            Nfine=5000
     93           Case(5)
    10094            file_id='/stellar_spectra/hd128167.txt'
    10195            tstellar=6700. ! Segura et al. (2003)
    102          else
     96            file_id_lam='/stellar_spectra/lam.txt'
     97            Nfine=5000
     98           Case(6)
     99            file_id='/stellar_spectra/BD_Teff-1600K.txt'
     100            tstellar=1600. ! Segura et al. (2003)
     101            file_id_lam='/stellar_spectra/lamBD.txt'
     102            Nfine=5000
     103           Case(7)
     104            file_id='/stellar_spectra/BD_Teff-1000K.txt'
     105            tstellar=1000. ! Segura et al. (2003)
     106            file_id_lam='/stellar_spectra/lamBD.txt'
     107            Nfine=5000
     108           Case Default
    103109            print*,'Error: unknown star type chosen'
    104110            call abort
     111         End Select
     112
     113         allocate(lam(Nfine),stel_f(Nfine))
     114
     115         file_path_lam=TRIM(datadir)//TRIM(file_id_lam)
     116         open(110,file=file_path_lam,form='formatted',status='old',iostat=ios)
     117         if (ios.ne.0) then        ! file not found
     118           write(*,*) 'Error from ave_stelspec'
     119           write(*,*) 'Data file ',trim(file_id_lam),' not found.'
     120           write(*,*)'Check that your path to datagcm:',trim(datadir)
     121           write(*,*)' is correct. You can change it in callphys.def with:'
     122           write(*,*)' datadir = /absolute/path/to/datagcm'
     123           write(*,*)' Also check that there is a ',trim(file_id_lam),' there.'
     124           call abort
     125         else
     126           do ifine=1,Nfine
     127             read(110,*) lam(ifine)
     128           enddo
     129           close(110)
    105130         endif
    106          
     131
     132
     133         ! load high resolution wavenumber data
    107134         file_path=TRIM(datadir)//TRIM(file_id)
    108135         open(111,file=file_path,form='formatted',status='old',iostat=ios)
    … …  
    121148           close(111)
    122149         endif
     150         
     151         ! sum data by band
     152         band=1
     153         Do while(lam(1).lt. real(10000.0/BWNV(band+1)))
     154            if (band.gt.L_NSPECTV) exit
     155            band=band+1
     156         enddo
     157         dl=lam(2)-lam(1)
     158         STELLAR(band)=STELLAR(band)+stel_f(1)*dl
     159         do ifine = 2,Nfine
     160            if(lam(ifine) .gt. real(10000.0/BWNV(band)))then
     161               band=band-1
     162            endif
     163            if(band .lt. 1) exit
     164            dl=lam(ifine)-lam(ifine-1)
     165            STELLAR(band)=STELLAR(band)+stel_f(ifine)*dl
     166         end do
     167               
     168         
     169         STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
     170         deallocate(lam,stel_f)
     171         
    123172      endif
    124173
    125 
    126       ! sum data by band
    127       band=1
    128       do ifine = 1,Nfine
    129 
    130          if(lam(Nfine-ifine+1) .lt. real(10000.0/BWNV(band+1)))then
    131             band=band+1
    132          endif
    133          if(band .gt. L_NSPECTV)then
    134             goto 9999 ! ok, ok, I know they're evil
    135          endif
    136          STELLAR(band)=STELLAR(band)+stel_f(Nfine-ifine+1)*dl
    137 
    138       end do
    139 
    140 9999  continue
    141       STELLAR=STELLAR/sum(STELLAR)
    142 
    143 
    144174      end subroutine ave_stelspec

  • trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90

    r650 r773  
    194194               hice(ig)        = max(hice(ig),0.0)
    195195               hice(ig)        = min(hice(ig),maxicethick)
    196                pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT* &
    197                    rhowater/pcapcal(ig)/ptimestep             
     196               pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep             
    198197               albedo(ig)      = albedoice
    199198

  • trunk/LMDZ.GENERIC/libf/phystd/largescale.F90

    r728 r773  
    9191      DO i = 1, ngridmx
    9292
    93          if(zt(i).le.15.) zt(i)=15.   ! check too low temperatures
     93         if(zt(i).le.15.) then
     94            print*,'in lsc',i,zt(i)
     95            zt(i)=15.   ! check too low temperatures
     96         endif
    9497!         call watersat(zt(i),pplay(i,k),zqs(i))
    9598         call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))

  • trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90

    r728 r773  
    5050!     Local variables
    5151      INTEGER i, k, iq, nn
    52       INTEGER, PARAMETER :: niter = 6
     52      INTEGER, PARAMETER :: niter = 1
    5353      INTEGER k1, k1p, k2, k2p
    5454      LOGICAL itest(ngridmx)

  • trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

    r751 r773  
    642642           if (tlocked) then
    643643              ztim1=SIN(declin)
    644               ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres)
    645               ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
     644!              ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres)
     645!              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
     646! JL12 corrects tidally resonant cases. nres=omega_rot/omega_orb
     647              ztim2=COS(declin)*COS(2.*pi*(zday/year_day)*(nres-1.))
     648              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day)*(nres-1.))
    646649
    647650              call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
    … …  
    10171020!             Moist convection
    10181021!             ----------------
    1019                ! Re-evaporate cloud water/ice
    1020 !              dqevap1(1:ngridmx,1:nlayermx)=0.
    1021 !              dtevap1(1:ngridmx,1:nlayermx)=0.
    1022 !               call evap(ptimestep,pt,pq,pdq,pdt,dqevap1,dtevap1,qevap,tevap)
    1023                
     1022
    10241023               dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0.
    10251024               dtmoist(1:ngridmx,1:nlayermx)=0.
    … …  
    10321031                          +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice)
    10331032               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx)
    1034 
    1035 !               do l=1,nlayermx
    1036 !                do ig=1,ngridmx
    1037 !                 if(isnan(dqmoist(ig,l,igcm_h2o_ice))) print*,'Nan in dqmoist ice, abort,ig,l',ig,l
    1038 !                 if(isnan(dqmoist(ig,l,igcm_h2o_vap))) print*,'Nan in dqmoist vap, abort,ig,l',ig,l
    1039 !                 if(isnan(dtmoist(ig,l))) print*,'Nan in dtmoist, abort,ig,l',ig,l               
    1040 !                  if(isnan(dqmoist(ig,l,igcm_h2o_ice)).or.isnan(dqmoist(ig,l,igcm_h2o_vap)).or.isnan(dtmoist(ig,l))) STOP
    1041 !                enddo
    1042 !              enddo
    1043 !              print*,'after moist',dqmoist(185,1,igcm_h2o_ice),dqmoist(185,1,igcm_h2o_vap),dtmoist(185,1)*86400.
    10441033
    10451034               !-------------------------
    … …  
    10511040                  enddo
    10521041                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
    1053                   print*,'In moistadj MAX atmospheric energy change   =',MAXVAL(pdt(:,:))*86400.,'K/day'
     1042                  print*,'In moistadj MAX atmospheric energy change   =',MAXVAL(dtmoist(:,:))*ptimestep,'K/step'
     1043                  print*,'In moistadj MIN atmospheric energy change   =',MINVAL(dtmoist(:,:))*ptimestep,'K/step'
    10541044                 
    10551045                ! test energy conservation

  • trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F

    r751 r773  
    299299!!! - physical constants: nevermind, things are done allright below
    300300!!! - physical frequency: nevermind, in inifis this is a simple print
     301      cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
    301302      CALL inifis(1,llm,day0,daysec,dtphys,
    302303     .            lati,long,area,rad,g,r,cpp)

  • trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90

    r716 r773  
    2424      use radcommon_h, only : tgasref,tgasmin,tgasmax
    2525      use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight
    26       use radcommon_h, only : wrefvar
     26      use radcommon_h, only : wrefvar,WNOI,WNOV
    2727      use datafile_mod, only: datadir
    2828
    … …  
    4949
    5050      real*8 x, xi(4), yi(4), ans, p
    51       real kappa_IR, kappa_VI
     51!     For gray case (JL12)
     52      real kappa_IR, kappa_VI, IR_VI_wnlimit, nVI_limit,nIR_limit
    5253
    5354      integer ngas, igas, jgas
    … …  
    274275!     Get gaseous k-coefficients and interpolate onto finer pressure grid
    275276
     277      if (graybody) then
     278!        constant absorption coefficient in visible
     279         write(*,*)"graybody: constant absorption coefficient in visible:"
     280         kappa_VI=-100000.
     281         call getin("kappa_VI",kappa_VI)
     282         write(*,*)" kappa_VI = ",kappa_VI
     283         kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27    ! conversion from m^2/kg to cm^2/molecule         
     284     
     285!        constant absorption coefficient in IR
     286         write(*,*)"graybody: constant absorption coefficient in InfraRed:"
     287         kappa_IR=-100000.
     288         call getin("kappa_IR",kappa_IR)
     289         write(*,*)" kappa_IR = ",kappa_IR       
     290         kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27    ! conversion from m^2/kg to cm^2/molecule
     291               
     292!        wavelength used to separate IR from VI
     293         IR_VI_wnlimit=3000.
     294         write(*,*)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
     295         
     296         nVI_limit=L_NSPECTV
     297         Do nw=1,L_NSPECTV
     298            if (WNOV(nw).gt.IR_VI_wnlimit) then
     299               nVI_limit=nw-1
     300               exit
     301            endif
     302         End do
     303         nIR_limit=L_NSPECTI
     304         Do nw=1,L_NSPECTI
     305            if (WNOI(nw).gt.IR_VI_wnlimit) then
     306               nIR_limit=nw-1
     307               exit
     308            endif
     309         End do
     310         write(*,*)"graybody: Visible / Infrared separation set at band ",nIR_limit,nVI_limit
     311             
     312      End if
     313
     314
    276315      ! VISIBLE
    277       if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not.graybody)) then
     316      if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then
    278317         file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
    279318         file_path=TRIM(datadir)//TRIM(file_id)
    … …  
    293332
    294333      else if (callgasvis.and.graybody) then
    295          ! constant absorption coefficient in visible
    296          write(*,*)"constant absorption coefficient in visible:"
    297          kappa_VI = -100000.
    298          call getin("kappa_VI",kappa_VI)
    299          write(*,*)" kappa_VI = ",kappa_VI
    300          kappa_VI = kappa_VI * 1.e4 * mugaz * 1.672621e-27       ! conversion from m^2/kg to cm^2/molecule         
    301          gasv8(:,:,:,:,:) = kappa_VI
    302 
     334         if(nVI_limit.eq.0) then
     335            gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_VI
     336         else if (nVI_limit.eq.L_NSPECTV) then
     337            gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_IR
     338         else
     339            gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)=kappa_IR
     340            gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)=kappa_VI
     341         end if
    303342      else
    304343         print*,'Visible corrk gaseous absorption is set to zero.'
    305          gasv8(:,:,:,:,:) = 0.0
    306 
     344         gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
    307345      endif
    308346
    309347      ! INFRA-RED
    310       if ((.not.TRIM(corrkdir).eq.'null').and.(.not.graybody)) then
     348      if ((.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then
    311349         file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
    312350         file_path=TRIM(datadir)//TRIM(file_id)
    … …  
    364402
    365403      else if (graybody) then
    366          ! constant absorption coefficient in IR
    367          write(*,*)"constant absorption coefficient in InfraRed:"
    368          kappa_IR = -100000.
    369          call getin("kappa_IR",kappa_IR)
    370          write(*,*)" kappa_IR = ",kappa_IR
    371          kappa_IR = kappa_IR * 1.e4 * mugaz * 1.672621e-27       ! conversion from m^2/kg to cm^2/molecule       
    372          gasi8(:,:,:,:,:) = kappa_IR
    373 
     404         if(nIR_limit.eq.0) then
     405            gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_VI
     406         else if (nIR_limit.eq.L_NSPECTI) then
     407            gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_IR
     408         else
     409            gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)=kappa_IR
     410            gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)=kappa_VI
     411         end if
    374412      else
    375413         print*,'Infrared corrk gaseous absorption is set to zero.'
    376          gasi8(:,:,:,:,:) = 0.0
    377 
     414         gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
    378415      endif
    379416
    … …  
    389426                  do nw=1,L_NSPECTV
    390427                     if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
    391                         gasv8(nt,np,nh,nw,ng) = &
    392                             log10(gasv8(nt,np,nh,nw,ng))
     428                        gasv8(nt,np,nh,nw,ng) = log10(gasv8(nt,np,nh,nw,ng))
    393429                     else
    394430                        gasv8(nt,np,nh,nw,ng) = -200.0
    … …  
    398434                  do nw=1,L_NSPECTI
    399435                     if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
    400                         gasi8(nt,np,nh,nw,ng) = &
    401                             log10(gasi8(nt,np,nh,nw,ng))
     436                        gasi8(nt,np,nh,nw,ng) = log10(gasi8(nt,np,nh,nw,ng))
    402437                     else
    403438                        gasi8(nt,np,nh,nw,ng) = -200.0
Note: See TracChangeset for help on using the changeset viewer.