Changeset 366


Ignore:
Timestamp:
Nov 9, 2011, 3:47:17 PM (13 years ago)
Author:
rwordsworth
Message:

OSR output bugs fixed.
Improvements to kcm for pure H2 atmospheres.

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r305 r366  
    3535#include "callkeys.h"
    3636#include "tracer.h"
     37#include "gases.h"
    3738
    3839!-----------------------------------------------------------------------
    … …  
    104105      REAL*8 tauaero(L_LEVELS+1,naerkind)
    105106      REAL*8 nfluxtopv,nfluxtopi,nfluxtop
    106       real*8 NFLUXTOPV_nu(L_NSPECTV)
    107       real*8 NFLUXTOPI_nu(L_NSPECTI)
     107      real*8 nfluxoutv_nu(L_NSPECTV) ! outgoing band-resolved VI flux at TOA (W/m2)
     108      real*8 nfluxtopi_nu(L_NSPECTI) ! net band-resolved IR flux at TOA (W/m2)
    108109      real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
    109110      REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
    … …  
    163164      logical OLRz
    164165      real OLR_nu(ngrid,L_NSPECTI)
    165       real GSR_nu(ngrid,L_NSPECTV)
     166      !real GSR_nu(ngrid,L_NSPECTV)
     167      real OSR_nu(ngrid,L_NSPECTV)
    166168      real*8 NFLUXGNDV_nu(L_NSPECTV)
    167169
    … …  
    184186      external CBRT
    185187
    186 !     included by RW for runway greenhouse 1D study
     188!     included by RW for runaway greenhouse 1D study
    187189      real muvar(ngridmx,nlayermx+1)
    188190      real vtmp(nlayermx)
    … …  
    589591      end do
    590592
     593
     594
     595!-----------------------------------------------------------------------
     596!     kcm mode only
    591597      if(kastprof)then
     598
     599         DO l=1,nlayer
     600            muvarrad(2*l)   = mugaz
     601            muvarrad(2*l+1) = mugaz
     602         END DO
     603
     604         do k=1,L_LEVELS
     605            qvar(k) = 0.0
     606         end do
     607         print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
     608      endif
     609
     610
     611      if(kastprof.and.(ngasmx.gt.1))then
    592612
    593613         DO l=1,nlayer
    … …  
    613633         qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O
    614634
    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
     635
     636      endif
     637
    621638
    622639      ! Keep values inside limits for which we have radiative transfer coefficients
    … …  
    720737
    721738            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
    722                  acosz,stel_fract,gweight,nfluxtopv,nfluxgndv_nu,  &
     739                 acosz,stel_fract,gweight,                         &
     740                 nfluxtopv,nfluxoutv_nu,nfluxgndv_nu,              &
     741                 !acosz,stel_fract,gweight,nfluxtopv,nfluxgndv_nu,  &
    723742                 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
    724743
    725744         else                          ! during the night, fluxes = 0
    726             nfluxtopv=0.0
     745            nfluxtopv       = 0.0
     746            nfluxoutv_nu(:) = 0.0
     747            nfluxgndv_nu(:) = 0.0
    727748            do l=1,L_NLAYRAD
    728749               fmnetv(l)=0.0
    … …  
    775796            end do
    776797            do nw=1,L_NSPECTV
    777                GSR_nu(ig,nw)=nfluxgndv_nu(nw)
     798               !GSR_nu(ig,nw)=nfluxgndv_nu(nw)
     799               OSR_nu(ig,nw)=nfluxoutv_nu(nw)
    778800            end do
    779801         endif
    … …  
    824846      if(specOLR)then
    825847        if(ngrid.ne.1)then
    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)
     848          !call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
     849          !call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
    828850        endif
    829851      endif
    … …  
    845867               close(117)
    846868
    847                open(127,file='GSRnu.out')
     869               open(127,file='OSRnu.out')
    848870               do nw=1,L_NSPECTV
    849                   write(127,*) GSR_nu(1,nw)
     871                  write(127,*) OSR_nu(1,nw)
    850872               enddo
    851873               close(127)

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

    r305 r366  
    4242     Pn = rho_n*T*rmn
    4343
    44      if(gnom(2).eq.'H2O')then
     44     if(ngasmx.eq.1)then
     45        print*,'Cannot have moist adiabat with one gas...'
     46        stop
     47     endif
     48
     49     if(gnom(ngasmx).eq.'H2O')then
    4550
    4651        call psat_H2O(T-2d-1,psat_minus)
    … …  
    6873        endif
    6974
    70      elseif(gnom(2).eq.'NH3')then
     75     elseif(gnom(ngasmx).eq.'NH3')then
    7176
    7277        call psat_NH3(T-2d-1,psat_minus)
    … …  
    104109  case(0) ! dry
    105110
    106      if(gnom(2).eq.'H2O')then
     111     cp_v=0.0
     112     if(gnom(ngasmx).eq.'H2O')then
    107113        cp_v = (32.24+1.923d-3*T+1.055d-5*T**2-3.511d-9*T**3)/m_v
    108      elseif(gnom(2).eq.'NH3')then
     114     elseif(gnom(ngasmx).eq.'NH3')then
    109115        cp_v = 2.058d3
    110      elseif(gnom(2).eq.'CH4')then
     116     elseif(gnom(ngasmx).eq.'CH4')then
    111117        cp_v = 2.226d3
    112118     endif

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

    r305 r366  
    162162  call su_gases
    163163  call calc_cpp_mugaz
     164
    164165
    165166  call inifis(1,llm,0,86400.0,1.0,0.0,0.0,1.0,rad,g,r,cpp)
    … …  
    192193           read(90,*,iostat=ierr) tnom(iq)
    193194           if (ierr.ne.0) then
    194               write(*,*) 'rcm1d: error reading tracer names...'
     195              write(*,*) 'kcm1d: error reading tracer names...'
    195196              stop
    196197           endif
    197198        enddo !of do iq=1,nq
    198199     endif
    199   endif
     200  !endif
    200201 
    201202  call initracer()
     203
     204  endif
     205
    202206
    203207  do iq=1,nqmx
    … …  
    215219
    216220  iter    = 1
    217   Tstrat  = 200.0
     221  Tstrat  = 60.0
    218222  dTstrat = 1000.0
    219223
    … …  
    227231     call kcmprof_fn(psurf,qsurf(1),tsurf,    &
    228232          tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
    229 
    230 
    231      !if(1.eq.2)then
    232233
    233234!    Run radiative transfer
    … …  
    242243     print*,'Tstrat = ',Tstrat
    243244     dTstrat = Tstrat
    244      Tstrat  = (fluxabs_sw/(2*sigma))**0.25 ! skin temperature (gray approx.)
     245     !Tstrat  = Tsurf*(0.2786*(psurf/100000.)**(-1.123))**0.25
     246     ! skin temperature (gray approx.) using analytic pure H2 expression
     247     !Tstrat  = (fluxabs_sw/(2*sigma))**0.25 ! skin temperature (gray approx.)
     248     Tstrat  = (fluxtop_lw/(2*sigma))**0.25 ! skin temperature (gray approx.)
    245249     dTstrat = dTstrat-Tstrat
    246 
    247      !endif
    248 
    249      !dTstrat = 0
    250250
    251251     if(abs(dTstrat).lt.1.0)then
    … …  
    264264
    265265! Calculate total atmospheric energy
    266   call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
    267      q(:,1),muvar,Eatmtot)
     266  Eatmtot=0.0
     267!  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
     268!     q(:,1),muvar,Eatmtot)
    268269
    269270! ------------------------

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

    r305 r366  
    4141  double precision Ptop, dlogp, Psat_max
    4242  parameter (Ptop=1.0)           ! Pressure at TOA [Pa]
    43   parameter (Psat_max=5000.0) ! Maximum vapour pressure [Pa]
     43  !parameter (Psat_max=100000.0) ! Maximum vapour pressure [Pa]
     44  parameter (Psat_max=0.0) ! set to zero for dry calculations
    4445
    4546  double precision T(1:nlay)                       ! temperature [K]
    … …  
    8384  ! modify/generalise later??
    8485
    85   if(gnom(2).eq.'H2O')then
     86  if(ngasmx.gt.2)then
     87     print*,'Only two species possible at present'
     88     stop
     89  elseif(ngasmx.eq.1)then
     90     if(psat_max.gt.0.0)then
     91        print*,'Must have Psat_max=0 if no variable species'
     92        stop
     93     endif
     94     print*, 'Assuming pure atmosphere'
     95     m_v   = 1.0
     96     tcrit = 1000.0
     97  elseif(gnom(ngasmx).eq.'H2O')then
    8698     m_v   = dble(mH2O/1000.)
    8799     tcrit = 6.47d2
    88   elseif(gnom(2).eq.'NH3')then
     100  elseif(gnom(ngasmx).eq.'NH3')then
    89101     m_v   = 17.031/1000.
    90102     tcrit = 4.06d2
    91   elseif(gnom(2).eq.'CH4')then
     103  elseif(gnom(ngasmx).eq.'CH4')then
    92104     m_v   = 16.04/1000.
    93105     tcrit = 1.91d2
    … …  
    113125!i!  endif
    114126
    115   if(gnom(2).eq.'H2O')then
     127
     128  psat_v=psat_max
     129  if(gnom(ngasmx).eq.'H2O')then
    116130     call Psat_H2O(tsurf,psat_v)
    117   elseif(gnom(2).eq.'NH3')then
     131  elseif(gnom(ngasmx).eq.'NH3')then
    118132     call Psat_NH3(tsurf,psat_v)
    119133  endif
    120134
    121   ! Moist adiabat unless greater than psat_max
     135  ! Moist adiabat unless greater than or equal to psat_max
    122136  if(psat_v*1d6.lt.psat_max)then
    123137    Psurf_v = Psat_v*1d6
    … …  
    141155  endif
    142156
     157
     158
    143159  ! define fine pressure grid
    144160  psurf_rcm = real(Psurf_n+Psurf_v)
    … …  
    211227
    212228     ! test for moist adiabat at next level
    213      if(gnom(2).eq.'H2O')then
     229     psat_v=psat_max
     230     if(gnom(ngasmx).eq.'H2O')then
    214231        call Psat_H2O(T(ilay+1),psat_v)
    215      elseif(gnom(2).eq.'NH3')then
     232     elseif(gnom(ngasmx).eq.'NH3')then
    216233        call Psat_NH3(T(ilay+1),psat_v)
    217234     endif
    … …  
    234251
    235252     if(profil_flag(ilay) .eq. 1)then
    236        
    237         if(gnom(2).eq.'H2O')then
     253       
     254        psat_v=psat_max
     255        if(gnom(ngasmx).eq.'H2O')then
    238256           call Psat_H2O(T(ilay+1),psat_v)
    239         elseif(gnom(2).eq.'NH3')then
     257        elseif(gnom(ngasmx).eq.'NH3')then
    240258           call Psat_NH3(T(ilay+1),psat_v)
    241259        endif

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

    r305 r366  
    157157            DCONT = DCONT*dz(k)
    158158
    159                if((.not.callgasvis))then
    160                   DCONT=0.0
    161                endif
     159            if((.not.callgasvis))then
     160               DCONT=0.0
     161            endif
    162162
    163163

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

    r253 r366  
    11      SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV,
    2      *                  UBAR0,STEL,GWEIGHT,NFLUXTOPV,NFLUXGNDV_nu,
     2     *                  UBAR0,STEL,GWEIGHT,NFLUXTOPV,NFLUXOUTV_nu,
     3     *                  NFLUXGNDV_nu,
    34     *                  FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf)
    45
    … …  
    1819      real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD)
    1920      real*8 NFLUXTOPV, FLUXUP, FLUXDN
    20       real*8 NFLUXTOPV_nu(L_NSPECTV)
     21      real*8 NFLUXOUTV_nu(L_NSPECTV)
    2122      real*8 NFLUXGNDV_nu(L_NSPECTV)
    2223      real*8 GWEIGHT(L_NGAUSS)
    … …  
    3839
    3940      DO NW=1,L_NSPECTV
    40          NFLUXTOPV_nu(NW)=0.0
     41         NFLUXOUTV_nu(NW)=0.0
    4142         NFLUXGNDV_nu(NW)=0.0
    4243      END DO
    … …  
    108109          END DO
    109110
    110 c     and same thing by spectral band... (RDW)
    111           NFLUXTOPV_nu(NW) = NFLUXTOPV_nu(NW)
    112      *      +(FLUXUP-FLUXDN)*GWEIGHT(NG)*
    113      *                          (1.0-FZEROV(NW))
     111c     band-resolved flux leaving TOA (RDW)
     112          NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
     113     *      +FLUXUP*GWEIGHT(NG)*(1.0-FZEROV(NW))
    114114
    115 
    116 c     flux at gnd (RDW)
     115c     band-resolved flux at ground (RDW)
    117116          NFLUXGNDV_nu(NW) = NFLUXGNDV_nu(NW)
    118117     *      +FMDV(L_NLAYRAD)*GWEIGHT(NG)*(1.0-FZEROV(NW))
    … …  
    163162        END DO
    164163
    165 c         and same thing by spectral band... (RDW)
    166           NFLUXTOPV_nu(NW) = NFLUXTOPV_nu(NW)
    167      *      +(FLUXUP-FLUXDN)*FZERO
     164c     band-resolved flux leaving TOA (RDW)
     165          NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
     166     *      +FLUXUP*FZERO
    168167
    169 c     flux at gnd (RDW)
     168c     band-resolved flux at ground (RDW)
    170169          NFLUXGNDV_nu(NW) = NFLUXGNDV_nu(NW)+FMDV(L_NLAYRAD)*FZERO
    171170
Note: See TracChangeset for help on using the changeset viewer.