Changeset 305


Ignore:
Timestamp:
Sep 22, 2011, 2:09:09 PM (13 years ago)
Author:
rwordsworth
Message:

Several new files added as part of the climate evolution model
(main program kcm.F90). Some general cleanup in physiq.F90 and
callcorrk.F90. Bugs in dust radiative transfer and H2 Rayleigh
scattering corrected.

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
14 added
5 deleted
14 edited

Legend:

Unmodified
Added
Removed

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

    r253 r305  
    181181!       2. Opacity calculation
    182182            if (aerofixed.or.(i_h2oice.eq.0).or.clearsky) then
    183                do ig=1, ngrid ! to stop the rad tran bug
     183               do ig=1,ngrid ! to stop the rad tran bug
    184184                  do l=1,nlayer
    185185                     aerosol(ig,l,iaer) =1.e-9
    186186                  end do
    187187               end do
     188
     189               ! put cloud at cloudlvl
     190               if(kastprof.and.(cloudlvl.ne.0.0))then
     191                  ig=1
     192                  do l=1,nlayer
     193                     if(int(cloudlvl).eq.l)then
     194                     !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
     195                        print*,'Inserting cloud at level ',l
     196                        !aerosol(ig,l,iaer)=10.0
     197
     198                        rho_ice=920.0
     199
     200                        ! the Kasting approximation
     201                        aerosol(ig,l,iaer) =                      &
     202                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
     203                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
     204                          !( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
     205                          ( 4.0e-4 + 1.E-9 ) *         &
     206                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
     207
     208
     209                        open(115,file='clouds.out',form='formatted')
     210                        write(115,*) l,aerosol(ig,l,iaer)
     211                        close(115)
     212
     213                        return
     214                     endif
     215                  end do
     216
     217                  call abort
     218               endif
     219
    188220            else
    189221               do ig=1, ngrid
    … …  
    232264!       2. Opacity calculation
    233265
    234             do l=1,nlayer
    235                do ig=1,ngrid-1 ! to stop the rad tran bug
     266            do ig=1,ngrid
     267               do l=1,nlayer
    236268                  zp=700./pplay(ig,l)
    237269                  aerosol(ig,l,1)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &

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

    r253 r305  
    4646
    4747      ! tau0/p0=tau/p (Hansen 1974)
    48       ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0)  *  4*delta^2/(g*mugaz*lambda^4)
    49       ! Then calculate tau0 = 1.16e-20 * 4*delta^2/(g*mugaz*lambda^4 [in um])
    50       ! where delta=n-1
     48      ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4)
     49      ! where delta=n-1 and N0 is an amagat
    5150
    5251      typeknown=.false.
    5352      do igas=1,ngasmx
    54          !if((igas.eq.vgas).or.(gfrac(igas).lt.1.e-4))then
    5553         if(igas.eq.vgas)then
    5654            print*,'Ignoring ',gnom(igas),' in Rayleigh scattering '// &
    … …  
    6058            'as its mixing ratio is less than 0.05.'
    6159            ! ignore variable gas in Rayleigh calculation
    62             ! ignore gases of mixing ratio < 1e-4 in Rayleigh calculation
     60            ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
    6361            tauconsti(igas) = 0.0
    6462         else
    … …  
    6765            elseif(gnom(igas).eq.'N2_')then
    6866               tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
     67            elseif(gnom(igas).eq.'H2O')then
     68               tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0
    6969            elseif(gnom(igas).eq.'H2_')then
    70                tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
     70               tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0
     71               !tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
    7172               ! uses n=1.000132 from Optics, Fourth Edition.
    72                ! around four times more opaque than CO2
    73                ! around 5.5 times more opaque than air
     73               ! was out by a factor 2!
    7474            elseif(gnom(igas).eq.'He_')then
    7575               print*,'Helium not ready yet!'
    … …  
    104104            tauvar=0.0
    105105            do igas=1,ngasmx
    106                !if((igas.eq.vgas).or.(gfrac(igas).lt.1.e-4))then
    107106               if((igas.eq.vgas).or.(gfrac(igas).lt.5.e-2))then
    108107                  ! ignore variable gas in Rayleigh calculation
    109                   ! ignore gases of mixing ratio < 1e-4 in Rayleigh calculation
     108                  ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
    110109                  tauvari(igas)   = 0.0
    111110               else
    … …  
    114113                  elseif(gnom(igas).eq.'N2_')then
    115114                     tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
     115                  elseif(gnom(igas).eq.'H2O')then
     116                     tauvari(igas) = 1.0/wl**4 ! to be improved...
    116117                  elseif(gnom(igas).eq.'H2_')then
    117                      tauvari(igas) = 1.0/wl**4 ! no wvl dependence of ref. index here - improve?
     118                     tauvari(igas) = 1.0/wl**4
    118119                  else
    119120                     print*,'Error in calc_rayleigh: Gas species not recognised!'
    … …  
    138139      end do
    139140
    140       print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-5), &
    141            ' um, tau_R = ',TAURAY(L_NSPECTV-5)*1013.25
    142       print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-6), &
    143            ' um, tau_R = ',TAURAY(L_NSPECTV-6)*1013.25
     141      print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-6),' um'
     142      print*,'tau_R = ',TAURAY(L_NSPECTV-6)*1013.25
     143      print*,'sig_R = ',TAURAY(L_NSPECTV-6)*g*mugaz*1.67e-27*100, &
     144             'cm^2 molecule^-1'
    144145
    145146    end subroutine calc_rayleigh

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

    r253 r305  
    1       subroutine callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,    &
     1      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,    &
    22          albedo,emis,mu0,pplev,pplay,pt,                      &
    3           tsurf,fract,dist_star,aerosol,cpp3D,                 &
     3          tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
    44          dtlw,dtsw,fluxsurf_lw,                               &
    55          fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,        &
    6           reffrad,tau_col,ptime,pday,cloudfrac,totcloudfrac,   &
     6          reffrad,tau_col,cloudfrac,totcloudfrac,   &
    77          clearsky,firstcall,lastcall)
     8
    89
    910      use radinc_h
    … …  
    117118      real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
    118119
    119       real*8 qvar(L_LEVELS)          ! mixing ratio of variable component
     120      real*8 qvar(L_LEVELS)          ! mixing ratio of variable component (mol/mol)
    120121      REAL pq(ngridmx,nlayer,nq)
    121122      REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
    … …  
    159160      character(len=2)  :: tempOLR
    160161      character(len=30) :: filenomOLR
    161       real ptime, pday
     162      !real ptime, pday
    162163      logical OLRz
    163164      real OLR_nu(ngrid,L_NSPECTI)
    164 
    165165      real GSR_nu(ngrid,L_NSPECTV)
    166166      real*8 NFLUXGNDV_nu(L_NSPECTV)
    … …  
    177177      real pqtest(ngridmx,nlayer,nq)
    178178
    179 !     for Dave Crisp LBL data comparison
    180       logical crisp
    181 
    182179!     are particle radii fixed?
    183180      logical radfixed
    … …  
    187184      external CBRT
    188185
    189       crisp=.false.
     186!     included by RW for runway greenhouse 1D study
     187      real muvar(ngridmx,nlayermx+1)
     188      real vtmp(nlayermx)
     189      REAL*8 muvarrad(L_LEVELS)
     190
    190191      radfixed=.false.
    191192
    … …  
    245246         Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
    246247
     248
    247249         if((igcm_h2o_vap.eq.0) .and. varactive)then
    248250            print*,'varactive in callcorrk but no h2o_vap tracer.'
    … …  
    265267      enddo
    266268
     269
     270      if(kastprof)then
     271         radfixed=.true.
     272      endif
    267273
    268274      if(radfixed)then
    … …  
    276282            do l=1,nlayer
    277283               do ig=1,ngrid
    278                   reffrad(ig,l,2) = 2.e-5 ! H2O ice
     284                  !reffrad(ig,l,2) = 2.e-5 ! H2O ice
     285                  reffrad(ig,l,2) = 5.e-6 ! H2O ice
    279286               enddo
    280287            enddo
    … …  
    297304
    298305         maxrad=0.0
     306         !averad=0.0
    299307         minrad=1.0
    300308         do l=1,nlayer
     309
     310            !masse = (pplev(ig,l) - pplev(ig,l+1))/g
     311
    301312            do ig=1,ngrid
    302 
    303313               if(tracer)then
    304314                  reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
    … …  
    308318               reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6)
    309319
     320               !averad = averad + reffrad(ig,l,1)*area(ig)
    310321               maxrad = max(reffrad(ig,l,1),maxrad)
    311322               minrad = min(reffrad(ig,l,1),minrad)
    … …  
    345356      endif
    346357
     358
    347359!     how much light we get
    348360      fluxtoplanet=0
    … …  
    360372           reffrad,QREFvis3d,QREFir3d,                             &
    361373           tau_col,cloudfrac,totcloudfrac,clearsky)                ! get aerosol optical depths
     374
    362375
    363376!-----------------------------------------------------------------------
    … …  
    512525      endif
    513526
    514       if(ngridmx.eq.1) then       ! fixed zenith angle in 1D
     527      if(ngridmx.eq.1) then       ! fixed zenith angle 'szangle' in 1D
    515528         acosz = cos(pi*szangle/180.0)
    516          !acosz=mu0(ig)
    517          !print*,'oned SZ angle OVERRIDE in callcorrk'
    518529         print*,'acosz=',acosz,', szangle=',szangle
    519530      else
    … …  
    521532      endif
    522533
    523 !-----------------------------------------------------------------------
    524 !     Water vapour (to be generalised / ignored for other planets)
     534
     535!-----------------------------------------------------------------------
     536!     Water vapour (to be generalised for other gases eventually)
    525537     
    526       if(varactive) then
     538      if(varactive)then
    527539
    528540         i_var=igcm_h2o_vap
    529 
    530541         do l=1,nlayer
    531 
    532542            qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
    533543            qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2   
    534544            ! Average approximation as for temperature...
    535 
    536545         end do
    537546         qvar(1)=qvar(2)
    … …  
    575584      end if
    576585
    577       ! Now convert from kg/kg to mol/mol
     586      ! IMPORTANT: Now convert from kg/kg to mol/mol
    578587      do k=1,L_LEVELS
    579588         qvar(k) = qvar(k)/epsi
    580589      end do
     590
     591      if(kastprof)then
     592
     593         DO l=1,nlayer
     594            muvarrad(2*l)   = muvar(ig,nlayer+2-l)
     595            muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
     596                                muvar(ig,max(nlayer+1-l,1)))/2
     597         END DO
     598     
     599         muvarrad(1) = muvarrad(2)
     600         muvarrad(2*nlayermx+1)=muvar(ig,1)
     601
     602         print*,'Recalculating qvar with VARIABLE epsi for kastprof'
     603         i_var=igcm_h2o_vap
     604         do l=1,nlayer
     605            vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O
     606         end do
     607
     608         do l=1,nlayer
     609            qvar(2*l)   = vtmp(nlayer+1-l)
     610            qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
     611         end do
     612         qvar(1)=qvar(2)
     613         qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O
     614
     615         !do k=1,L_LEVELS
     616         !   qvar(k) = 1.0
     617         !end do
     618         !print*,'ASSUMING qH2O=1 EVERYWHERE IN CALLCORRK!!'
     619
     620      endif
    581621
    582622      ! Keep values inside limits for which we have radiative transfer coefficients
    … …  
    607647      tlevrad(2*nlayermx+1)=tsurf(ig)
    608648     
    609 
    610       if(crisp)then
    611          open(111,file='/u/rwlmd/CrispLBL/p.dat')
    612          open(112,file='/u/rwlmd/CrispLBL/T.dat')
    613          do k=1,L_LEVELS
    614             read(111,*) plevrad(k)
    615             read(112,*) tlevrad(k)
    616             plevrad(k)=plevrad(k)*1000.0
    617          enddo
    618          close(111)
    619          close(112)
    620       endif
    621 
    622649      tmid(1) = tlevrad(2)
    623650      tmid(2) = tlevrad(2)
    … …  
    651678            print*,'Minimum temperature is outside the radiative'
    652679            print*,'transfer kmatrix bounds, exiting.'
    653             !print*,'WARNING, OVERRIDING FOR TEST'
    654             call abort
     680            print*,'WARNING, OVERRIDING FOR TEST'
     681            !call abort
    655682         elseif(tlevrad(k).gt.tgasmax)then
    656683            print*,'Maximum temperature is outside the radiative'
    657684            print*,'transfer kmatrix bounds, exiting.'
    658             !print*,'WARNING, OVERRIDING FOR TEST'
    659             call abort
     685            print*,'WARNING, OVERRIDING FOR TEST'
     686            !print*, 'T=',pt
     687            !call abort
    660688         endif
    661689      enddo
    … …  
    688716            call optcv(dtauv,tauv,taucumv,plevrad,                 &
    689717                 qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
    690                  tmid,pmid,taugsurf,qvar)
     718                 tmid,pmid,taugsurf,qvar,muvarrad)
    691719
    692720
    … …  
    704732         end if
    705733
    706         !print*,'taugsurf=',taugsurf(:,L_NGAUSS-1)
     734
    707735
    708736
    … …  
    712740         call optci(plevrad,tlevrad,dtaui,taucumi,                  &
    713741              qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
    714               taugsurfi,qvar)
     742              taugsurfi,qvar,muvarrad)
    715743
    716744         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
    … …  
    729757         fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
    730758         fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
     759
    731760
    732761         if(fluxtop_dn(ig).lt.0.0)then
    … …  
    740769         endif
    741770
    742          if(crisp)then
    743             open(111,file='/u/rwlmd/CrispLBL/GCMfluxdn.dat')
    744             open(112,file='/u/rwlmd/CrispLBL/GCMfluxup.dat')
    745             do k=1,L_NLAYRAD
    746                write(111,*) fluxdni(k)
    747                write(112,*) fluxupi(k)
    748             enddo
    749             close(111)
    750             close(112)
    751             print*,'fluxdni',fluxdni
    752             print*,'fluxupi',fluxupi
    753             call abort
    754          endif
    755 
    756771!     Spectral output, for exoplanet observational comparison
    757772         if(specOLR)then
    … …  
    809824      if(specOLR)then
    810825        if(ngrid.ne.1)then
    811           call writediagspec(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
    812 !          call writediagspec(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
     826          call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
     827          call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W m^-2",3,GSR_nu)
    813828        endif
    814829      endif
    815830
    816       if(ngrid.eq.1)then
    817          open(113,file='surf_vals_long.out',ACCESS='append')
    818          write(113,*) tsurf(1),fluxtop_dn(1),         &
    819               real(-nfluxtopv),real(nfluxtopi)
    820          close(113)
    821       endif
    822 
    823       !if(kastprof)then ! for kasthop only
    824 
    825         if(ngrid.eq.1)then
    826 
    827            print*,'Saving tsurf,psurf in surf_vals.out...'
     831      if(lastcall.and.(ngrid.eq.1))then
     832
     833           print*,'Saving scalar quantities in surf_vals.out...'
     834           print*,'psurf = ', pplev(1,1),' Pa'
    828835           open(116,file='surf_vals.out')
    829836           write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1),         &
    … …  
    843850               enddo
    844851               close(127)
    845            endif
    846 
    847 !     mixing ratio: CO2 ice particles
    848            if(igcm_co2_ice.ne.0)then
    849               open(122,file='qCO2ice.out')
    850               do l=1,L_NLAYRAD
    851                  write(122,*) pq(1,l,igcm_co2_ice)
    852               enddo
    853               close(122)
    854            endif
    855 
    856 !     mixing ratio: H2O vapour
    857            if(igcm_h2o_vap.ne.0)then
    858               open(123,file='qH2Ovap.out')
    859               do l=1,L_NLAYRAD
    860                  write(123,*) pq(1,l,igcm_h2o_vap)
    861               enddo
    862               close(123)
    863852           endif
    864853
    … …  
    878867              close(119)
    879868           endif
    880         endif
    881 
    882       !endif
     869
     870      endif
    883871
    884872    end subroutine callcorrk

  • trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

    r253 r305  
    1616     &   , kastprof                                                     &
    1717     &   , noradsurf                                                    &
    18      &   , Tsurfk                                                       &
     18     &   , graybody                                                     &
    1919     &   , Tstrat                                                       &
    2020     &   , newtonian                                                    &
    … …  
    4545     &   , n2mixratio                                                   &
    4646     &   , co2supsat                                                    &
     47     &   , cloudlvl                                                     &
    4748     &   , pceil                                                   
    4849
    … …  
    7374      logical CLFvarying
    7475      logical noradsurf
     76      logical graybody
    7577
    7678      integer iddist
    … …  
    9799      real maxicethick
    98100      real Tsaldiff
    99       real Tsurfk
    100101      real tau_relax
     102      real cloudlvl

  • trunk/LMDZ.GENERIC/libf/phystd/gases.h

    r253 r305  
    88      real gfrac(ngasmx)
    99
    10       COMMON/gases/gnom,gfrac,vgas
     10      COMMON/gases/gnom,vgas,gfrac
    1111
    1212!-----------------------------------------------------------------------

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

    r253 r305  
    5353      parameter (snowlayer=33.0)        ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
    5454      real oceantime
    55       parameter (oceantime=1*3600)
     55      parameter (oceantime=10*24*3600)
    5656
    5757      logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
    … …  
    101101
    102102      oceanbulkavg=.false.
    103       activerunoff=.false.
     103      activerunoff=.true.
    104104      oceanalbvary=.false.
    105105
    … …  
    311311         END DO
    312312
    313          print*,'Mean ocean temperature = ',tsea
     313         print*,'Mean ocean temperature               = ',tsea,' K'
    314314
    315315      endif

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

    r253 r305  
    1313!
    1414!   Initialisation for the physical parametrisations of the LMD
    15 !   martian atmospheric general circulation modele.
     15!   Generic Model.
    1616!
    1717!   author: Frederic Hourdin 15 / 10 /93
    … …  
    5656
    5757
     58
    5859      REAL prad,pg,pr,pcpp,pdaysec,ptimestep
    5960 
    … …  
    265266         write(*,*)" kastprof = ",kastprof
    266267
    267          write(*,*)"Surface temperature in kastprof mode?"
    268          Tsurfk=273.15
    269          call getin("Tsurfk",Tsurfk)
    270          write(*,*)" Tsurfk = ",Tsurfk
     268         write(*,*)"Uniform absorption coefficient in IR?"
     269         graybody=.false.
     270         call getin("graybody",graybody)
     271         write(*,*)" graybody = ",graybody
    271272
    272273! Test of incompatibility:
    … …  
    344345         write(*,*)" Fat1AU = ",Fat1AU
    345346
    346          write(*,*)"Set temperature to 1 K above CO2 condensation?"
    347          nearco2cond=.false.
    348          call getin("nearco2cond",nearco2cond)
    349          write(*,*)" nearco2cond = ",nearco2cond
    350347
    351348!     1D solar zenith angle
    … …  
    394391         endif
    395392
     393         write(*,*)"Cloud pressure level (with kastprof only):"
     394         cloudlvl=0. ! default value
     395         call getin("cloudlvl",cloudlvl)
     396         write(*,*)" cloudlvl = ",cloudlvl
     397
    396398         write(*,*)"Is the variable gas species radiatively active?"
     399         Tstrat=167.0
    397400         varactive=.false.
    398401         call getin("varactive",varactive)

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

    r253 r305  
    2727      integer nS,nT
    2828      parameter(nS=1649)
    29 !      parameter(nT=7)
    3029      parameter(nT=14)
    3130
    … …  
    3433      double precision temp_arr(nT)
    3534      double precision abs_arr(nS,nT)
    36       !double precision data_tmp(nT+1)
    3735      double precision data_tmp(nT/2+1)
    3836
    … …  
    5250
    5351         ! cold
    54          dt_file=datafile(1:LEN_TRIM(datafile))//'/H2H2_CIA_LT.dat'
     52         dt_file=datafile(1:LEN_TRIM(datafile))//'/continuum_data/H2H2_CIA_LT.dat'
    5553         open(33,file=dt_file,form='formatted',status='old',iostat=ios)
    5654         if (ios.ne.0) then        ! file not found
    … …  
    6967
    7068         ! hot
    71          dt_file=datafile(1:LEN_TRIM(datafile))//'/H2H2_CIA_HT.dat'
     69         dt_file=datafile(1:LEN_TRIM(datafile))//'/continuum_data/H2H2_CIA_HT.dat'
    7270         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
    7371         if (ios.ne.0) then        ! file not found
    … …  
    106104         print*,'   pressure ',pres,' Pa'
    107105
    108          call bilinear(wn_arr,temp_arr,nS,nT,abs_arr,wn,temp,abcoef)
     106         call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
    109107
    110108         print*,'the absorption is ',abcoef,' cm^-1 amg^-2'
    … …  
    117115      else
    118116 
    119          call bilinear(wn_arr,temp_arr,nS,nT,abs_arr,wn,temp,abcoef)
     117         call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
    120118         abcoef=abcoef*100.0*amagat**2 ! convert to m^-1
    121119         !print*,'The absorption is ',abcoef,' m^-1'
    … …  
    130128
    131129
    132 
    133 ! put as separate file eventually
    134 
    135130!-------------------------------------------------------------------------
    136       subroutine bilinear(x_arr,y_arr,nX,nY,f2d_arr,x_in,y_in,f)
     131      subroutine bilinearH2H2(x_arr,y_arr,f2d_arr,x_in,y_in,f)
    137132!     Necessary for interpolation of continuum data
    138133
    … …  
    140135
    141136      integer nX,nY,i,j,a,b
     137      parameter(nX=1649)
     138      parameter(nY=14)
    142139
    143140      real*8 x_in,y_in,x,y,x1,x2,y1,y2
    … …  
    176173      if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
    177174         write(*,*) 'Warning from bilinear.for:'
    178          write(*,*) 'Extrapolating outside H2-H2 CIA temperature range!'
     175         write(*,*) 'Outside continuum temperature range!'
    179176         if(y.lt.y_arr(1))then
    180177            y=y_arr(1)+0.01
    … …  
    212209     
    213210      return
    214     end subroutine bilinear
     211    end subroutine bilinearH2H2

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

    r253 r305  
    11      subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
    22           QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
    3            TMID,PMID,TAUGSURF,QVAR)
     3           TMID,PMID,TAUGSURF,QVAR,MUVAR)
    44
    55      use radinc_h
    … …  
    5959      real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
    6060      real*8 DCONT
    61       double precision wn_cont, p_cont, T_cont
    62 
    63 !     Water mixing ratio variables
    64       real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS)
     61      double precision wn_cont, p_cont, p_air, T_cont, dtemp
     62
     63!     Variable species mixing ratio variables
     64      real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
    6565      real*8  KCOEF(4)
    6666      integer NVAR(L_LEVELS)
    … …  
    9797      do K=2,L_LEVELS
    9898         DPR(k) = PLEV(K)-PLEV(K-1)
    99          U(k)   = Cmk*DPR(k)    ! only Cmk line in optci.F
    10099
    101100         !--- Kasting's CIA ----------------------------------------
    … …  
    110109
    111110         ! if we have continuum opacities, we need dz
    112          dz(k)=dpr(k)*R*TMID(K)/(g*PMID(K))
     111         if(kastprof)then
     112            dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
     113            U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k)
     114         else
     115            dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
     116            U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
     117         endif
    113118
    114119         call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
    … …  
    132137            DCONT = 0.0 ! continuum absorption
    133138
    134             ! include H2 continuum
     139            ! include continua if necessary
     140            wn_cont = dble(wnoi(nw))
     141            T_cont  = dble(TMID(k))
    135142            do igas=1,ngasmx
    136                if(gnom(igas).eq.'H2_')then
    137 
    138                   wn_cont = dble(wnoi(nw))
    139                   T_cont  = dble(TMID(k))
    140                   p_cont  = dble(PMID(k)*scalep*gfrac(igas))
    141 
    142                   call interpolateH2H2(wn_cont,T_cont,p_cont,DCONT,.false.)
    143 
    144                   DCONT=DCONT*dz(k)
    145 
    146 
     143
     144               if(gfrac(igas).eq.-1)then ! variable
     145                  p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
     146               else
     147                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
    147148               endif
     149
     150               dtemp=0.0
     151               if(gnom(igas).eq.'H2_')then
     152                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
     153               elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
     154                  p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
     155                  call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
     156
     157               endif
     158
     159               DCONT = DCONT + dtemp
     160
    148161            enddo
     162
     163
     164            DCONT = DCONT*dz(k)
    149165
    150166            !--- Kasting's CIA ----------------------------------------
    … …  
    154170            !----------------------------------------------------------
    155171
    156             ! Water continuum currently inactive!
    157             !if(varactive)then
    158             !   call water_cont(PMID(K),TMID(K),WRATIO(K),WNOI(nw),DCO2)
    159             !endif
    160172
    161173            do ng=1,L_NGAUSS-1
    … …  
    197209               
    198210               TAUGAS  = U(k)*ANS
     211
     212               if(graybody)then
     213                  TAUGAS = 0.0
     214                  DCONT  = U(k)*3.3e-26
     215               endif
    199216
    200217               TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
    … …  
    311328      END DO
    312329
    313 
    314330      return
     331
     332
    315333    end subroutine optci
    316334
    317335
    318 !-------------------------------------------------------------------------
    319     subroutine water_cont(p,T,wratio,nu,DCONT)
    320 !   Calculates the continuum opacity for H2O
    321 !   NOT CURRENTLY USED
    322 
    323       implicit none
    324 
    325       real p, T, wratio, nu, DCONT
    326       real h1, h2
    327 
    328       h1 = exp(1800.*(1./T - 0.0034))
    329       h2 = 1.25e-22 + 1.67e-19*exp(-2.62e-13*nu)
    330 
    331 !      DCONT = h1*h2*p*(p*wratio)**2/(R*T)
    332 !      DCONT=0.0 ! to be implemented...
    333 
    334       return
    335 
    336     end subroutine water_cont
    337 
     336

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

    r253 r305  
    11      SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
    22          QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
    3           TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR)
     3          TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
    44
    55      use radinc_h
    … …  
    3030!     TLEV(L) - Temperature at the layer boundary
    3131!     PLEV(L) - Pressure at the layer boundary (i.e. level)
    32 !     GASV(NT,NPS,NW,NG) - Visual CO2 k-coefficients
     32!     GASV(NT,NPS,NW,NG) - Visible k-coefficients
    3333!     
    3434!-------------------------------------------------------------------
    … …  
    6565      real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER
    6666
    67 !     Water mixing ratio variables
    68       real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS)
     67!     Variable species mixing ratio variables
     68      real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
    6969      real*8 KCOEF(4)
    7070      integer NVAR(L_LEVELS)
    … …  
    7474
    7575!     variables for k in units m^-1
    76       double precision wn_cont, p_cont, T_cont
     76      double precision wn_cont, p_cont, p_air, T_cont, dtemp
    7777      real*8 dz(L_LEVELS), DCONT
    7878
    … …  
    9292         DPR(k) = PLEV(K)-PLEV(K-1)
    9393
    94          U(k)   = Cmk*DPR(k)
    95 
    9694         ! if we have continuum opacities, we need dz
    97          dz(k)  = dpr(k)*R*TMID(K)/(g*PMID(K))
     95         if(kastprof)then
     96            dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
     97            U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k)
     98         else
     99            dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
     100            U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
     101         endif
    98102
    99103         call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
    … …  
    128132            DCONT = 0.0 ! continuum absorption
    129133
    130             ! include H2 continuum
     134            ! include continua if necessary
     135            wn_cont = dble(wnov(nw))
     136            T_cont  = dble(TMID(k))
    131137            do igas=1,ngasmx
     138
     139               if(gfrac(igas).eq.-1)then ! variable
     140                  p_cont  = dble(PMID(k)*scalep*QVAR(K)) ! qvar = mol/mol
     141               else
     142                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
     143               endif
     144
     145               dtemp=0.0
    132146               if(gnom(igas).eq.'H2_')then
    133 
    134                   wn_cont = dble(wnov(nw))
    135                   T_cont  = dble(TMID(k))
    136                   p_cont  = dble(PMID(k)*scalep*gfrac(igas))
    137                   call interpolateH2H2(wn_cont,T_cont,p_cont,DCONT,.false.)
    138 
    139                   DCONT=DCONT*dz(k)
    140 
    141                   if((.not.callgasvis))then
    142                      DCONT=0.0
    143                   endif
    144 
    145 
    146                endif
     147                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
     148               elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
     149                  p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
     150                  call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
     151               endif
     152
     153               DCONT = DCONT + dtemp
     154
    147155            enddo
     156
     157            DCONT = DCONT*dz(k)
     158
     159               if((.not.callgasvis))then
     160                  DCONT=0.0
     161               endif
     162
    148163
    149164            do NG=1,L_NGAUSS-1
    … …  
    184199
    185200               TAUGAS          = U(k)*ANS
     201
     202               if(graybody)then
     203                  TAUGAS = 0.0
     204                  DCONT  = 0.0
     205               endif
     206
     207
    186208
    187209               !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS

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

    r253 r305  
    534534      IF (ierr.NE.NF_NOERR) THEN
    535535         PRINT*, 'phyetat0: Lecture echouee pour <totcloudfrac>'
    536          CALL abort
     536         !CALL abort
    537537      ENDIF
    538538      xmin = 1.0E+20
    … …  
    541541      xmax = MAXVAL(totcloudfrac)
    542542      PRINT*,'Cloud fraction <totcloudfrac>:', xmin, xmax
    543 
    544543
    545544

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

    r253 r305  
    88 
    99      use radinc_h, only : naerkind
    10       use watercommon_h, only : mx_eau_sol, To, RLVTT
     10      use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O
    1111     
    1212      implicit none
    … …  
    392392      integer igtest
    393393
    394 !     included by RW for pressure broadening broadening study
    395       integer gascut
     394!     included by RW for runway greenhouse 1D study
     395      real muvar(ngridmx,nlayermx+1)
    396396
    397397!     included by RW for variable H2O particle sizes
    … …  
    404404!     included by RW for sourceevol
    405405      real ice_initial(ngridmx)!, isoil
     406      real delta_ice,ice_tot
    406407      real ice_min(ngridmx)
    407408      save ice_min
    … …  
    502503            close(128)
    503504       
    504             if(num_run.ne.0.and.mod(num_run,6).eq.0)then
     505            if(num_run.ne.0.and.mod(num_run,3).eq.0)then
    505506               print*,'Updating ice at end of this year!'
    506507               ice_update=.true.
    … …  
    713714
    714715              if(kastprof)then
    715                  ISR=fract(1)*Fat1AU/dist_star**2/2.0
    716 
    717                  call simpleprof_fn(tsurf,pt,zpopsk,ISR,pplev(1,1))
    718                  open(116,file='profpres.out',form='formatted')
    719                  open(117,file='proftemp.out',form='formatted')
    720                  do l=1,nlayer
    721                     write(116,*) pplay(1,l)
    722                     write(117,*) pt(1,l)
    723                  enddo
    724                  close(116)
    725                  close(117)
     716                 print*,'kastprof should not = true here'
     717                 call abort
    726718              endif
    727 
     719              muvar(:,:)=0.0 ! only used for climate evolution for now
    728720
    729721!             Option to call scheme once for clear regions, once for cloudy regions (BC)
    … …  
    731723
    732724                 clearsky=.true.
    733                  call callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,           &
     725                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
    734726                      albedo,emis,mu0,pplev,pplay,pt,                      &
    735                       tsurf,fract,dist_star,aerosol,cpp3D,                 &
     727                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
    736728                      zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
    737                       fluxabs_sw1,fluxtop_dn,reffrad,tau_col1,ptime,pday,  &
     729                      fluxabs_sw1,fluxtop_dn,reffrad,tau_col1,             &
    738730                      cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
    739731
    740732                 clearsky=.false.
    741                  call callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,           &
     733                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
    742734                      albedo,emis,mu0,pplev,pplay,pt,                      &
    743                       tsurf,fract,dist_star,aerosol,cpp3D,                 &
     735                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
    744736                      zdtlw2,zdtsw2,fluxsurf_lw2,fluxsurf_sw2,fluxtop_lw2, &
    745                       fluxabs_sw2,fluxtop_dn,reffrad,tau_col2,ptime,pday,  &
     737                      fluxabs_sw2,fluxtop_dn,reffrad,tau_col2,             &
    746738                      cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
    747739
    … …  
    767759
    768760                 clearsky=.false.
    769                  call callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,         &
     761                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,         &
    770762                      albedo,emis,mu0,pplev,pplay,pt,                    &
    771                       tsurf,fract,dist_star,aerosol,cpp3D,               &
     763                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,         &
    772764                      zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
    773                       fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday,  &
     765                      fluxabs_sw,fluxtop_dn,reffrad,tau_col,             &
    774766                      cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
    775767
    … …  
    881873      endif ! of if (callrad)
    882874
    883        
    884875!-----------------------------------------------------------------------
    885876!    4. Vertical diffusion (turbulent mixing):
    … …  
    906897              zdum1,zdum2,zdh,pdq,zflubid,           &
    907898              zdudif,zdvdif,zdhdif,zdtsdif,q2,       &
    908               zdqdif,zdqsdif)
     899              zdqdif,zdqsdif,lastcall)
    909900
    910901         do l=1,nlayer
    … …  
    11431134         endif
    11441135
    1145          call condens_co2cloud(ngrid,nlayer,nq,ptimestep,   &
     1136         call condense_cloud(ngrid,nlayer,nq,ptimestep,   &
    11461137              capcal,pplay,pplev,tsurf,pt,                  &
    11471138              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
    … …  
    12131204               call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
    12141205               DO l = 1, nlayer 
    1215                   DO ig = 1, ngrid              
     1206                  DO ig = 1, ngrid
    12161207                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
    12171208                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
    … …  
    16261617
    16271618!     Increment surface temperature
    1628       if(.not.kastprof)then
    1629          do ig=1,ngrid
    1630             tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
    1631          enddo
    1632       endif
     1619      do ig=1,ngrid
     1620         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
     1621      enddo
    16331622
    16341623!     Compute soil temperatures and subsurface heat flux
    … …  
    17831772         endif
    17841773
    1785          if(kastprof)then ! for f(p,T) computation (not needed in universal model)
    1786             open(15,file='ene_bal.out',form='formatted',access='append')
    1787             write(15,*) pplev(1,1), Tsurf(1), ISR, OLR, ASR
    1788             close(15)
    1789             open(16,file='Seff.out',form='formatted')
    1790             write(16,*) OLR/ASR
    1791             close(16)
    1792             open(17,file='alb.out',form='formatted')
    1793             write(17,*) 1.0-ASR/ISR
    1794             close(17)
    1795             call abort
    1796          endif
    1797 
    17981774      endif
    17991775
    … …  
    18301806         do ig=1,ngrid
    18311807            do l=1,nlayer
    1832                reffcol(ig,1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * reffrad(ig,l,1) *    &
    1833                     (pplev(ig,l) - pplev(ig,l+1)) / g
     1808               if(co2cond)then
     1809                  reffcol(ig,1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * &
     1810                                  reffrad(ig,l,1) *    &
     1811                                  (pplev(ig,l) - pplev(ig,l+1)) / g
     1812               endif
    18341813               if(water)then
    1835                reffcol(ig,2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * reffrad(ig,l,2) *    &
    1836                     (pplev(ig,l) - pplev(ig,l+1)) / g
    1837 
     1814                  reffcol(ig,2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * &
     1815                                  reffrad(ig,l,2) *    &
     1816                                  (pplev(ig,l) - pplev(ig,l+1)) / g
    18381817               endif
    18391818            enddo
    … …  
    19231902!           Update surface ice distribution to iterate to steady state if requested
    19241903            if(ice_update)then
     1904
    19251905               do ig = 1, ngrid
    19261906
    1927                   ! if ice amount is declining, set it to zero
    1928                   if(qsurf_hist(ig,igcm_h2o_ice).lt.ice_initial(ig))then
    1929                      qsurf_hist(ig,igcm_h2o_ice)=0.0
    1930                      qsurf_hist(ig,igcm_h2o_vap)=0.0
     1907                  delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
     1908
     1909                  ! add one hundred years/timesteps of evolution
     1910                  qsurf_hist(ig,igcm_h2o_ice) = &
     1911                     qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0
     1912
     1913                  ! if ice has gone -ve, set to zero
     1914                  if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
     1915                     qsurf_hist(ig,igcm_h2o_ice) = 0.0
     1916                     qsurf_hist(ig,igcm_h2o_vap) = 0.0
    19311917                  endif
    1932  
    1933                   ! if ice amount is increasing and was never less than 0.5 kg/m2,
    1934                   ! set value in gridbox to 1x10^3 kg/m2
    1935                   !if(qsurf_hist(ig,igcm_h2o_ice).gt.0.0 .and. ice_min(ig).gt.0.5)then
    1936                   if(qsurf_hist(ig,igcm_h2o_ice).gt.ice_initial(ig))then
    1937                      if(ice_min(ig).gt.0.5)then ! ignore seasonal ice
    1938                         qsurf_hist(ig,igcm_h2o_ice)=1.e3
    1939                         qsurf_hist(ig,igcm_h2o_vap)=0.0
    1940                      endif
     1918
     1919                  ! set maximum to ice
     1920                  if(qsurf_hist(ig,igcm_h2o_ice).gt.1.e3)then
     1921                     qsurf_hist(ig,igcm_h2o_ice) = 1.e3
     1922                     qsurf_hist(ig,igcm_h2o_vap) = 0.0
    19411923                  endif
    19421924
    19431925               enddo
     1926
     1927               ! enforce ice conservation
     1928               ice_tot=0.0
     1929               do ig = 1, ngrid
     1930                  ice_tot = ice_tot + qsurf_hist(ig,igcm_h2o_ice)*area(ig)
     1931               enddo
     1932               do ig = 1, ngrid
     1933                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice)*(icesrf/ice_tot)
     1934               enddo
     1935
     1936
     1937
    19441938            endif
    19451939

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

    r253 r305  
    7474
    7575!      integer, parameter :: L_NPREF   = 8
     76!      integer, parameter :: L_NTREF   = 11
     77!      integer, parameter :: L_REFVAR  = 8   ! N2_H2Ovar
     78
     79!      integer, parameter :: L_NPREF   = 8
    7680!      integer, parameter :: L_NTREF   = 5
    7781!      integer, parameter :: L_REFVAR  = 1   ! therm_test2
    … …  
    109113      ! equivalent temperatures are 1/10 of these values
    110114      integer, parameter :: NTstar = 500
    111 !      integer, parameter :: NTstop = 9000 ! for hotstuff runs
    112       integer, parameter :: NTstop = 6000 ! for GJ581d / earlymars runs
    113 !      integer, parameter :: NTstop = 10000 ! for H2 warming runs
     115      integer, parameter :: NTstop = 9000 ! new default for all non hot Jupiter runs
     116      !integer, parameter :: NTstop = 6000 ! for GJ581d / earlymars runs
    114117
    115118! Maximum number of grain size classes

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

    r253 r305  
    493493
    494494      if(autozlevs)then
    495          if(kastprof)then
    496             Hscale=10.0
    497             Hmax=400.0
    498          else
    499495           
    500             open(91,file="z2sig.def",form='formatted')
    501             read(91,*) Hscale
    502             DO ilayer=1,nlayer-2
    503                read(91,*) Hmax
    504             enddo
    505             close(91)
    506  
    507          endif
     496         open(91,file="z2sig.def",form='formatted')
     497         read(91,*) Hscale
     498         DO ilayer=1,nlayer-2
     499            read(91,*) Hmax
     500         enddo
     501         close(91)
     502 
    508503           
    509504         print*,'Hmax = ',Hmax,' km'
Note: See TracChangeset for help on using the changeset viewer.