Ignore:
Timestamp:
Feb 13, 2012, 6:30:26 PM (13 years ago)
Author:
jleconte
Message:

13/02/2012 == JL + AS

  • All outputs are now in netCDF format. Even in 1D (No more G1D)
  • Clean up of the call to callcorrk when CLFvarying=true
  • Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
  • Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
  • Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
  • corrects a small bug when no clouds are present in aeropacity
Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
17 edited

Legend:

Unmodified
Added
Removed

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

    r486 r526  
    235235                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
    236236                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
    237                           ( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
    238                           ( pplev(ig,l) - pplev(ig,l+1) ) / g /   &
    239                             CLF1
     237                            pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
     238                          ( pplev(ig,l) - pplev(ig,l+1) ) / g /   & !   opacity in the clearsky=true and the
     239                            CLF1                                    !   clear=false/pq=0 case
    240240
    241241                     aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0))

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

    r486 r526  
    1       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,    &
     1      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
    22          albedo,emis,mu0,pplev,pplay,pt,                      &
    33          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,cloudfrac,totcloudfrac,   &
     6          OLR_nu,OSR_nu,                                       &
     7          reffrad,tau_col,cloudfrac,totcloudfrac,              &
    78          clearsky,firstcall,lastcall)
    89
    … …  
    8384      REAL fluxabs_sw(ngridmx)    ! SW flux absorbed by planet (W/m2)
    8485      REAL fluxtop_dn(ngridmx)    ! incident top of atmosphere SW flux (W/m2)
    85 
     86      REAL OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
     87      REAL OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
    8688!-----------------------------------------------------------------------
    8789!     Declaration of the variables required by correlated-k subroutines
    … …  
    154156!      real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
    155157
    156 !     for OLR spec
    157       integer OLRcount
    158       save OLRcount
    159       integer OLRcount2
    160       save OLRcount2
    161       character(len=2)  :: tempOLR
    162       character(len=30) :: filenomOLR
    163158      !real ptime, pday
    164159      logical OLRz
    165       real OLR_nu(ngrid,L_NSPECTI)
    166       !real GSR_nu(ngrid,L_NSPECTV)
    167       real OSR_nu(ngrid,L_NSPECTV)
    168160      real*8 NFLUXGNDV_nu(L_NSPECTV)
    169161
    … …  
    272264            stop
    273265         endif
     266         
     267         OLR_nu=0.
     268         OSR_nu=0.
    274269
    275270         firstcall=.false.   
    … …  
    331326
    332327            do ig=1,ngrid
    333                if(tracer)then
     328               if(tracer.and.igcm_co2_ice.gt.0)then
    334329                  reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
    335330                       (4*Nmix_co2*pi*rho_co2) )
    … …  
    561556         end do
    562557         qvar(1)=qvar(2)
    563          qvar(2*nlayermx+1)=qsurf(ig,i_var)
     558!         qvar(2*nlayermx+1)=qsurf(ig,i_var) !JL12 not very good to compare kg/kg and kg/m2???
    564559
    565560      elseif(varfixed)then
    … …  
    590585         Ttemp = tsurf(ig)
    591586         call watersat(Ttemp,ptemp,qsat)
     587         
    592588
    593589         !qvar(2*nlayermx+1)=qsat      ! fully saturated everywhere
    594590         qvar(2*nlayermx+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
     591
     592
     593!!!!!!!!!!!!!!!!!!!!!!!!  JL: for completely constant water profile uncoment the following line
     594!        qvar=0.005
    595595
    596596      else
    … …  
    766766              qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
    767767              taugsurfi,qvar,muvarrad)
    768 
     768             
    769769         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
    770770              wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, &
    … …  
    797797         if(specOLR)then
    798798            do nw=1,L_NSPECTI
    799                OLR_nu(ig,nw)=nfluxtopi_nu(nw)
     799               OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
    800800            end do
    801801            do nw=1,L_NSPECTV
    802802               !GSR_nu(ig,nw)=nfluxgndv_nu(nw)
    803                OSR_nu(ig,nw)=nfluxoutv_nu(nw)
     803               OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
    804804            end do
    805805         endif
    … …  
    848848
    849849!     IR spectral output, for exoplanet observational comparison
    850       if(specOLR)then
    851         if(ngrid.ne.1)then
    852           call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
    853           call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
    854         endif
    855       endif
    856 
    857       if(lastcall.and.(ngrid.eq.1))then
     850
     851
     852      if(lastcall.and.(ngrid.eq.1))then  ! could disable the 1D output, they are in the diagfi and diagspec... JL12
    858853
    859854           print*,'Saving scalar quantities in surf_vals.out...'
    … …  
    864859           close(116)
    865860
    866            if(specOLR)then
    867                open(117,file='OLRnu.out')
    868                do nw=1,L_NSPECTI
    869                   write(117,*) OLR_nu(1,nw)
    870                enddo
    871                close(117)
    872 
    873                open(127,file='OSRnu.out')
    874                do nw=1,L_NSPECTV
    875                   write(127,*) OSR_nu(1,nw)
    876                enddo
    877                close(127)
    878            endif
     861!           if(specOLR)then
     862!               open(117,file='OLRnu.out')
     863!               do nw=1,L_NSPECTI
     864!                  write(117,*) OLR_nu(1,nw)
     865!               enddo
     866!               close(117)
     867!
     868!               open(127,file='OSRnu.out')
     869!               do nw=1,L_NSPECTV
     870!                  write(127,*) OSR_nu(1,nw)
     871!               enddo
     872!               close(127)
     873!           endif
    879874
    880875!     OLR vs altitude: do it as a .txt file

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

    r486 r526  
    99     &   , enertest                                                     &
    1010     &   , callgasvis                                                   &
     11     &   , Continuum                                                    &
    1112     &   , Nmix_co2, Nmix_h2o                                           &
    1213     &   , dusttau                                                      &
    … …  
    5253     &   , season,diurnal,tlocked,lwrite                                &
    5354     &   , callstats,calleofdump                                        &
    54      &   , callgasvis
     55     &   , callgasvis,Continuum
    5556
    5657      logical enertest

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

    r486 r526  
    4646!     Francois Forget (1996)
    4747!     Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
     48!     Includes simplifed nucleation by J. Leconte (2011)
    4849!     
    4950!==================================================================
    … …  
    9697      REAL zcpi
    9798      REAL ztcond (ngridmx,nlayermx)
     99      REAL ztnuc (ngridmx,nlayermx)
    98100      REAL ztcondsol(ngridmx)
    99101      REAL zdiceco2(ngridmx)
    … …  
    271273!      end if
    272274
    273 !     Forecast the atmospheric frost temperature 'ztcond'
     275!     Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'
    274276      DO l=1,nlayer
    275277         DO ig=1,ngrid
    276278            ppco2=gfrac(igasco2)*pplay(ig,l)
    277279            call get_tcond_co2(ppco2,ztcond(ig,l))
     280            call get_tnuc_co2(ppco2,ztnuc(ig,l))
    278281         ENDDO
    279282      ENDDO
    … …  
    379382
    380383
    381                IF ((zt(ig,l).LT.ztcond(ig,l)).or.(zq(ig,l,i_co2ice).gt.0)) THEN
     384   !     ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011)
     385               IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_co2ice).gt.1.E-10)) THEN
    382386                  condsub(ig)=.true.
    383387                  pdtc(ig,l)   = (ztcond(ig,l) - zt(ig,l))/subptimestep
    … …  
    524528      real, parameter :: ptriple=518000.0
    525529
    526       peff=p!/co2supsat
     530      peff=p
    527531
    528532      if(peff.lt.ptriple)then
    … …  
    535539
    536540    end subroutine get_tcond_co2
     541   
     542   
     543   
     544 
     545!-------------------------------------------------------------------------
     546    subroutine get_tnuc_co2(p,tnuc)
     547!   Calculates the nucleation temperature for CO2, based on a simple super saturation criterion
     548!     (JL 2011)
     549
     550      implicit none
     551
     552#include "callkeys.h"
     553
     554      real p, peff, tnuc
     555      real, parameter :: ptriple=518000.0
     556
     557      peff=p/co2supsat
     558
     559      if(peff.lt.ptriple)then
     560         tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
     561      else
     562         tnuc = 684.2-92.3*log(peff)+4.32*log(peff)**2
     563         ! liquid-vapour transition (based on CRC handbook 2003 data)
     564      endif
     565      return
     566
     567    end subroutine get_tnuc_co2

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

    r486 r526  
    192192         call getin("callgasvis",callgasvis)
    193193         write(*,*) " callgasvis = ",callgasvis
     194
     195         write(*,*) "call continuum opacities in radiative transfer ?",
     196     &              "(matters only if callrad=T)"
     197         Continuum=.false. ! default value
     198         call getin("Continuum",Continuum)
     199         write(*,*) " Continuum = ",Continuum
    194200         
    195201         write(*,*) "call turbulent vertical diffusion ?"

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

    r305 r526  
    22
    33      use radinc_h, only: L_NSPECTI
    4       use radcommon_h, only: WAVEI
     4      use radcommon_h, only: WNOI
    55
    66      implicit none
    … …  
    5454      integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
    5555 !     integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
    56       integer :: idim_bandsIR ! "bandsIR" dimension ID #
     56      integer :: idim_bandsIR ! "IR Wavenumber" dimension ID #
    5757      integer, dimension(2) :: id 
    5858
    … …  
    114114      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
    115115!      ierr = NF_DEF_DIM (nid, "subsurface_layers",nsoilmx,idim_nsoilmx)
    116       ierr = NF_DEF_DIM (nid, "bandsIR",L_NSPECTI,idim_bandsIR)
     116      ierr = NF_DEF_DIM (nid, "IR Wavenumber",L_NSPECTI,idim_bandsIR)
    117117
    118118      ierr = NF_ENDDEF(nid)
    … …  
    290290      ! define variable
    291291#ifdef NC_DOUBLE
    292       ierr=NF_DEF_VAR(nid,"bandsIR",NF_DOUBLE,1,idim_bandsIR,nvarid)
    293 #else
    294       ierr=NF_DEF_VAR(nid,"bandsIR",NF_FLOAT,1,idim_bandsIR,nvarid)
    295 #endif
    296       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 27,
    297      .                        "Band limits in the infrared")
    298       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"cm^-1")
     292      ierr=NF_DEF_VAR(nid,"IR Wavenumber",NF_DOUBLE,1,
     293     .                          idim_bandsIR,nvarid)
     294#else
     295      ierr=NF_DEF_VAR(nid,"IR Wavenumber",NF_FLOAT,1,
     296     .                          idim_bandsIR,nvarid)
     297#endif
     298      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 34,
     299     .                        "Band mid frequency in the infrared")
     300      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",5,"cm^-1")
    299301      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
    300302      ! write variable
    301303#ifdef NC_DOUBLE
    302       ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WAVEI))
    303 #else
    304       ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WAVEI))
     304      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WNOI))
     305#else
     306      ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WNOI))
    305307#endif
    306308

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

    r305 r526  
    22
    33      use radinc_h, only: L_NSPECTV
    4       use radcommon_h, only: WAVEV
     4      use radcommon_h, only: WNOV
    55
    66      implicit none
    … …  
    5454      integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
    5555 !     integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
    56       integer :: idim_bandsVI ! "bandsIR" dimension ID #
     56      integer :: idim_bandsVI ! "VI Wavenumber" dimension ID #
    5757      integer, dimension(2) :: id 
    5858
    … …  
    114114      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
    115115!      ierr = NF_DEF_DIM (nid, "subsurface_layers",nsoilmx,idim_nsoilmx)
    116       ierr = NF_DEF_DIM (nid, "bandsVI",L_NSPECTV,idim_bandsVI)
     116      ierr = NF_DEF_DIM (nid, "VI Wavenumber",L_NSPECTV,idim_bandsVI)
    117117
    118118      ierr = NF_ENDDEF(nid)
    … …  
    290290      ! define variable
    291291#ifdef NC_DOUBLE
    292       ierr=NF_DEF_VAR(nid,"bandsVI",NF_DOUBLE,1,idim_bandsVI,nvarid)
    293 #else
    294       ierr=NF_DEF_VAR(nid,"bandsVI",NF_FLOAT,1,idim_bandsVI,nvarid)
    295 #endif
    296       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 27,
    297      .                        "Band limits in the visible")
    298       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"cm^-1")
     292      ierr=NF_DEF_VAR(nid,"VI Wavenumber",NF_DOUBLE,1,
     293     .                          idim_bandsVI,nvarid)
     294#else
     295      ierr=NF_DEF_VAR(nid,"VI Wavenumber",NF_FLOAT,1,
     296     .                          idim_bandsVI,nvarid)
     297#endif
     298      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 33,
     299     .                        "Band mid frequency in the visible")
     300      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",5,"cm^-1")
    299301      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
    300302      ! write variable
    301303#ifdef NC_DOUBLE
    302       ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WAVEV))
    303 #else
    304       ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WAVEV))
     304      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WNOV))
     305#else
     306      ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WNOV))
    305307#endif
    306308

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

    r374 r526  
    150150
    151151      else
    152  
     152
    153153         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
    154154         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
    … …  
    162162         abcoef = abcoef*Nmolec/(100.0**2) ! convert to m^-1
    163163         abcoef = abcoef*(presS/(presF+presS))     ! take H2O mixing ratio into account
     164
    164165
    165166         ! unlike for Rayleigh scattering, we do not currently weight by the BB function

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

    r471 r526  
    163163
    164164            DCONT = DCONT*dz(k)
     165           
     166            if(.not.Continuum)then
     167               DCONT=0.0
     168            endif
    165169
    166170            !--- Kasting's CIA ----------------------------------------
    … …  
    217221               TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
    218222               !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
    219 
    220223               DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
    221224

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

    r471 r526  
    157157            DCONT = DCONT*dz(k)
    158158
    159             if((.not.callgasvis))then
     159            if(.not.(callgasvis.and.Continuum))then
    160160               DCONT=0.0
    161161            endif
    162 
    163162
    164163            do NG=1,L_NGAUSS-1

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

    r486 r526  
    77                  pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
    88 
    9       use radinc_h, only : naerkind
     9      use radinc_h, only : naerkind,L_NSPECTI,L_NSPECTV
    1010      use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O
    1111      use gases_h
    … …  
    116116#include "comsaison.h"
    117117#include "control.h"
    118 #include "comg1d.h"
    119118#include "tracer.h"
    120119#include "watercap.h"
    … …  
    190189
    191190      character*80 fichier
    192       integer l,ig,ierr,iq,i, tapphys
     191      integer l,ig,ierr,iq,i, tapphys,nw
    193192
    194193      real fluxsurf_lw(ngridmx)      ! incident LW (IR) surface flux (W.m-2)
    … …  
    199198      real fluxtop_dn(ngridmx)
    200199      real fluxdyn(ngridmx)          ! horizontal heat transport by dynamics     
    201       save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn
     200      real OLR_nu(ngridmx,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
     201      real OSR_nu(ngridmx,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
     202      save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu
    202203
    203204
    … …  
    326327      real zpopskNI(ngridmx,nlayermx)
    327328
    328 !     included by RW to make 1D saves not every timestep
    329       integer countG1D
    330       save countG1D
    331329
    332330!     included by BC for evaporation     
    … …  
    357355!     included by BC for double radiative transfer call
    358356      logical clearsky
    359 
    360       real zdtsw1(ngridmx,nlayermx),zdtsw2(ngridmx,nlayermx)
    361       real zdtlw1(ngridmx,nlayermx),zdtlw2(ngridmx,nlayermx)
    362       real fluxsurf_lw1(ngridmx), fluxsurf_lw2(ngridmx)     
    363       real fluxsurf_sw1(ngridmx), fluxsurf_sw2(ngridmx)   
    364       real fluxtop_lw1(ngridmx), fluxtop_lw2(ngridmx)     
    365       real fluxabs_sw1(ngridmx), fluxabs_sw2(ngridmx)   
    366       real tau_col1(ngrid), tau_col2(ngrid)
     357      real zdtsw1(ngridmx,nlayermx)
     358      real zdtlw1(ngridmx,nlayermx)
     359      real fluxsurf_lw1(ngridmx)     
     360      real fluxsurf_sw1(ngridmx) 
     361      real fluxtop_lw1(ngridmx)   
     362      real fluxabs_sw1(ngridmx) 
     363      real tau_col1(ngrid)
     364      real OLR_nu1(ngrid,L_NSPECTI)
     365      real OSR_nu1(ngrid,L_NSPECTV)
    367366      real tf, ntf
    368367
    … …  
    423422      if (firstcall) then
    424423
    425          if(ngrid.eq.1)then
    426             saveG1D=day_step*10
    427             !saveG1D=1
    428             !saveG1D=-1 ! to never save
    429             countG1D=1
    430          endif
    431424
    432425!        variables set to 0
    … …  
    447440!        read startfi (initial parameters)
    448441!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    449 
     442          print*,'in physiqu i am in phyetat0'
    450443         call phyetat0("startfi.nc",0,0,nsoilmx,nq,           &
    451444               day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
    … …  
    681674
    682675      if (callrad) then
    683          if( mod(icount-1,iradia).eq.0) then
     676         if( mod(icount-1,iradia).eq.0.or.lastcall) then
    684677
    685678!          Compute local stellar zenith angles
    … …  
    720713              endif
    721714              muvar(:,:)=0.0 ! only used for climate evolution for now
    722 
    723 !             Option to call scheme once for clear regions, once for cloudy regions (BC)
     715             
     716             
     717!             standard callcorrk
     718              clearsky=.false.
     719              call callcorrk(ngrid,nlayer,pq,nq,qsurf,         &
     720                      albedo,emis,mu0,pplev,pplay,pt,                    &
     721                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,         &
     722                      zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
     723                      fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
     724                      reffrad,tau_col,cloudfrac,totcloudfrac,            &
     725                      clearsky,firstcall,lastcall)     
     726             
     727
     728!             Option to call scheme once more for clear regions
    724729              if(CLFvarying)then
    725730
    726                !!! THIS IS NOT A VERY GOOD STRATEGY HERE (everything should be in callcorrk)
    727                !!! ---> PROBLEMS WITH ALLOCATED ARRAYS
    728                !!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
    729 
     731                 !!! ---> PROBLEMS WITH ALLOCATED ARRAYS
     732                 !!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
    730733                 clearsky=.true.
    731734                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
    … …  
    733736                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
    734737                      zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
    735                       fluxabs_sw1,fluxtop_dn,reffrad,tau_col1,             &
    736                       cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
    737 
    738                  clearsky=.false.
    739                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
    740                       albedo,emis,mu0,pplev,pplay,pt,                      &
    741                       tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
    742                       zdtlw2,zdtsw2,fluxsurf_lw2,fluxsurf_sw2,fluxtop_lw2, &
    743                       fluxabs_sw2,fluxtop_dn,reffrad,tau_col2,             &
    744                       cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
     738                      fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
     739                      reffrad,tau_col1,cloudfrac,totcloudfrac,             &
     740                      clearsky,firstcall,lastcall)
     741                 clearsky = .false.  !! just in case.     
    745742
    746743                 ! Sum the fluxes and heating rates from cloudy/clear cases
    … …  
    749746                    ntf=1.-tf
    750747                   
    751                     fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw2(ig)
    752                     fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw2(ig)
    753                     fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw2(ig)
    754                     fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw2(ig)
    755                     tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col2(ig)
     748                    fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
     749                    fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
     750                    fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
     751                    fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
     752                    tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
    756753                   
    757754                    do l=1,nlayer
    758                        zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw2(ig,l)
    759                        zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw2(ig,l)
     755                       zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw(ig,l)
     756                       zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw(ig,l)
    760757                    enddo
     758
     759                    do nw=1,L_NSPECTV
     760                       OSR_nu(ig,nw) = ntf*OSR_nu1(ig,nw) + tf*OSR_nu(ig,nw)                   
     761                    enddo
     762                    do nw=1,L_NSPECTI
     763                       OLR_nu(ig,nw) = ntf*OLR_nu1(ig,nw) + tf*OLR_nu(ig,nw)                   
     764                    enddo
     765
    761766                 enddo
    762767
    763 !             Otherwise standard callcorrk
    764               else
    765 
    766                  clearsky=.false.
    767                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,         &
    768                       albedo,emis,mu0,pplev,pplay,pt,                    &
    769                       tsurf,fract,dist_star,aerosol,cpp3D,muvar,         &
    770                       zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
    771                       fluxabs_sw,fluxtop_dn,reffrad,tau_col,             &
    772                       cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
    773 
    774               endif
     768              endif !CLFvarying
    775769
    776770!             Radiative flux from the sky absorbed by the surface (W.m-2)
    … …  
    872866            print*,'In corrk SW atmospheric energy change =',dEtotSW,' W m-2'
    873867            print*,'In corrk LW atmospheric energy change =',dEtotLW,' W m-2'
     868            print*,'atmospheric energy change   (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
    874869            print*,'In corrk SW surface energy change     =',dEtotsSW,' W m-2'
    875870            print*,'In corrk LW surface energy change     =',dEtotsLW,' W m-2'
     871            print*,'surface energy change       (SW+LW)   =',dEtotsLW+dEtotsSW,' W m-2'
    876872         endif
    877873!-------------------------
    … …  
    884880
    885881      if (calldifv) then
    886 
     882     
    887883         do ig=1,ngrid
    888884            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
    … …  
    17431739            OLR = OLR + area(ig)*fluxtop_lw(ig)
    17441740            GND = GND + area(ig)*fluxgrd(ig)
     1741            if(.not.callsoil) GND=GND+ area(ig)*fluxrad(ig)
    17451742            DYN = DYN + area(ig)*fluxdyn(ig)
    17461743           
    … …  
    17651762         if(enertest)then
    17661763            print*,'SW energy balance SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
     1764            print*,'LW energy balance LW++ + ***ASR*** = ',dEtotLW+dEtotsLW+ASR/totarea,' W m-2'
    17671765         endif
    17681766
    17691767         if(meanOLR)then
    1770             if((ngridmx.gt.1) .or. (countG1D.eq.saveG1D))then
     1768            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
    17711769               ! to record global radiative balance
    17721770               open(92,file="rad_bal.out",form='formatted',access='append')
    … …  
    18571855
    18581856         if(meanOLR)then
    1859             if((ngridmx.gt.1) .or. (countG1D.eq.saveG1D))then
     1857            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
    18601858               ! to record global water balance
    18611859               open(98,file="h2o_bal.out",form='formatted',access='append')
    … …  
    18911889
    18921890
    1893       if (ngrid.ne.1) then
    18941891         print*,''
    18951892         print*,'--> Ls =',zls*180./pi
    … …  
    19771974!                        "Solar radiative flux to space","W.m-2",2,         &
    19781975!                        fluxtop_sw_tot)
     1976
     1977            call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
     1978            call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
     1979            call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
     1980
    19791981            call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
    19801982            call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
    … …  
    19841986
    19851987           if (tracer) then
     1988             do iq=1,nq
     1989                call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
     1990                call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
     1991                          'kg m^-2',2,qsurf(1,iq) )
     1992
     1993                call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
     1994                          'kg m^-2',2,qcol(1,iq) )
     1995                call wstats(ngridmx,trim(noms(iq))//'_reff',                          &
     1996                          trim(noms(iq))//'_reff',                                   &
     1997                          'm',3,reffrad(1,1,iq))
     1998              end do
    19861999            if (water) then
    19872000               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
    … …  
    20412054
    20422055!     Output aerosols
    2043         call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))
    2044         call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))
    2045         call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))
    2046         call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))
     2056        if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))
     2057        if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))
     2058        if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))
     2059        if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))
    20472060
    20482061!     Output tracers
    … …  
    20932106       endif
    20942107
    2095 !        ----------------------------------------------------------
    2096 !        Output in netcdf file "diagsoil.nc" for subterranean
    2097 !        variables (output every "ecritphy", as for writediagfi)
    2098 !        ----------------------------------------------------------
    2099 
    2100        ! Write soil temperature
    2101        !call writediagsoil(ngrid,"soiltemp","Soil temperature","K",3,tsoil)
    2102 
    2103       else         ! if(ngrid.eq.1)
    2104 
    2105 !        ----------------------------------------------------------------------
    2106 !        Output in grads file "g1d" (ONLY when using testphys1d)
    2107 !        ----------------------------------------------------------------------
    2108 
    2109          if(countG1D.eq.saveG1D)then
    2110          
    2111             call writeg1d(ngrid,1,tsurf,'tsurf','K')
    2112             call writeg1d(ngrid,1,ps,'ps','Pa')
    2113             call writeg1d(ngrid,nlayer,zt,'T','K')
    2114             call writeg1d(ngrid,nlayer,pq(1,1,1),'q','kg/kg')
    2115             call writeg1d(ngrid,nlayer,aerosol,'aerosol','SI')
    2116             call writeg1d(ngrid,nlayer,reffrad,'reffrad','SI')
    2117             call writeg1d(ngrid,nlayer,zdtlw,'dtlw','SI')
    2118             call writeg1d(ngrid,nlayer,zdtsw,'dtsw','SI')
    2119             call writeg1d(ngrid,nlayer,zdtdyn,'dtdyn','SI')
    2120 
    2121 !           radiation balance
    2122             call writeg1d(ngrid,1,ISR/totarea,'ISR','W m-2')
    2123             call writeg1d(ngrid,1,ASR/totarea,'ASR','W m-2')
    2124             call writeg1d(ngrid,1,OLR/totarea,'OLR','W m-2') 
    2125 
    2126             if(tracer) then
    2127                do iq=1,nq
    2128                   call writeg1d(ngrid,1,qsurf(1,iq),trim(noms(iq))//'_s','kg m^-2')
    2129                   call writeg1d(ngrid,1,qcol(1,iq),trim(noms(iq))//'_c','kg m^-2')
    2130                   call writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
    2131                end do
    2132                
    2133                if(waterrain)then
    2134                   CALL writeg1d(ngrid,1,zdqsrain,'rainfall','kg m-2 s-1')
    2135                   CALL writeg1d(ngrid,1,zdqssnow,'snowfall','kg m-2 s-1')
    2136                endif
    2137                
    2138             end if
    2139 
    2140             countG1D=1
    2141          else
    2142             countG1D=countG1D+1
    2143          endif ! if time to save
    2144 
    2145       endif        ! if(ngrid.ne.1)
     2108!      output spectrum
     2109      if(specOLR.and.corrk)then     
     2110                 call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
     2111                 call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
     2112      endif
     2113
    21462114
    21472115      icount=icount+1

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

    r253 r526  
    6262
    6363      logical simple                    ! Use very simple Emanuel scheme?
    64       parameter(simple=.true.)
     64      parameter(simple=.false.)
    6565
    6666      logical evap_prec                 ! Does the rain evaporate?

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

    r486 r526  
    4141#include "comvert.h"
    4242#include "netcdf.inc"
    43 #include "comg1d.h"
    4443#include "watercap.h"
    4544#include "fisice.h"
    … …  
    117116c=======================================================================
    118117
    119       saveprofile=.true.
     118      saveprofile=.false.
    120119
    121120c ------------------------------------------------------
    … …  
    274273      write(*,*)" time = ",time
    275274      time=time/24.E+0 ! convert time (hours) to fraction of sol
     275     
     276       ecritphy=day_step ! default value for ecritphy
     277       PRINT *,'Nunber of steps between writediagfi ?'
     278       call getin("ecritphy",ecritphy)
     279       write(*,*) " ecritphy = ",ecritphy
     280
    276281
    277282c  Discretisation (Definition de la grille et des pas de temps)
    … …  
    294299      ndt=ndt*day_step     
    295300      dtphys=daysec/day_step 
     301
     302
     303c output spectrum?
     304      write(*,*)"Output spectral OLR?"
     305      specOLR=.false.
     306      call getin("specOLR",specOLR)
     307      write(*,*)" specOLR = ",specOLR
    296308
    297309
    … …  
    304316      write(*,*) " psurf = ",psurf
    305317c Pression de reference
    306       pa=20.
    307       preff=610. ! these values are not needed in 1D 
     318      pa=psurf/30.
     319      preff=psurf ! these values are not needed in 1D  (are you sure JL12)
    308320
    309321c Gravity of planet
    … …  
    577589      enddo
    578590
    579 c    Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
    580 c    ----------------------------------------------------------------
    581 c    (effectuee avec "writeg1d" a partir de "physiq.F"
    582 c    ou tout autre subroutine physique, y compris celle ci).
    583 
    584         g1d_nlayer=nlayer
    585         g1d_nomfich='g1d.dat'
    586         g1d_unitfich=40
    587         g1d_nomctl='g1d.ctl'
    588         g1d_unitctl=41
    589         g1d_premier=.true.
    590         g2d_premier=.true.
    591591
    592592c  Ecriture de "startfi"
    … …  
    613613         lastcall=.true.
    614614         call stellarlong(day*1.0,zls)
    615          write(103,*) 'Ls=',zls*180./pi
    616          write(103,*) 'Lat=', lati(1)*180./pi
    617          write(103,*) 'Tau=', tauvis/700*psurf
    618          write(103,*) 'RunEnd - Atmos. Temp. File'
    619          write(103,*) 'RunEnd - Atmos. Temp. File'
    620          write(104,*) 'Ls=',zls*180./pi
    621          write(104,*) 'Lat=', lati(1)
    622          write(104,*) 'Tau=', tauvis/700*psurf
    623          write(104,*) 'RunEnd - Atmos. Temp. File'
     615!         write(103,*) 'Ls=',zls*180./pi
     616!         write(103,*) 'Lat=', lati(1)*180./pi
     617!         write(103,*) 'Tau=', tauvis/700*psurf
     618!         write(103,*) 'RunEnd - Atmos. Temp. File'
     619!         write(103,*) 'RunEnd - Atmos. Temp. File'
     620!         write(104,*) 'Ls=',zls*180./pi
     621!         write(104,*) 'Lat=', lati(1)
     622!         write(104,*) 'Tau=', tauvis/700*psurf
     623!         write(104,*) 'RunEnd - Atmos. Temp. File'
    624624        ENDIF
    625625
    … …  
    748748c    ========================================================
    749749
    750 c    fermeture pour conclure les sorties format grads dans "g1d.dat"
    751 c    et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
    752 c    ou tout autre subroutine physique, y compris celle ci).
    753 
    754750! save temperature profile
    755751      if(saveprofile)then
    … …  
    761757      endif
    762758
    763       logplevs(1)=log(plev(1)/100000.0)
    764       do ilayer=2,nlayer
    765          logplevs(ilayer)=log(play(ilayer)/100000.0)
    766       enddo
    767       CALL endg1d(1,nlayer,logplevs,ndt) ! now simply log(p) with p in bars
    768759
    769760c    ========================================================

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

    r486 r526  
    6464C     TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
    6565
    66       TTOP  = TLEV(2)
     66      TTOP  = TLEV(2)  ! JL12 why not (1) ???
    6767      TSURF = TLEV(L_LEVELS)
    6868

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

    r323 r526  
    196196        ! are undefined; so set them to 1
    197197        iphysiq=1
    198         irythme=1
     198        !JL11 irythme=1
    199199        ! NB:
    200200      endif

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

    r305 r526  
    112112! The following test is here to enforce that writediagfi is not used with the
    113113! 1D version of the GCM
    114       if (ngrid.eq.1) return
     114!not anymore (JL12)
     115      if (ngrid.eq.-1) return
    115116     
    116117c     nom=trim((nom))
    … …  
    238239              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
    239240              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
    240               ierr= NF_INQ_DIMID(nid,"bandsIR",id(3))
     241              ierr= NF_INQ_DIMID(nid,"IR Wavenumber",id(3))
    241242              ierr= NF_INQ_DIMID(nid,"Time",id(4))
    242243

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

    r305 r526  
    112112! The following test is here to enforce that writediagfi is not used with the
    113113! 1D version of the GCM
    114       if (ngrid.eq.1) return
     114!not anymore (JL12)
     115      if (ngrid.eq.-1) return
    115116     
    116117c     nom=trim((nom))
    … …  
    238239              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
    239240              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
    240               ierr= NF_INQ_DIMID(nid,"bandsVI",id(3))
     241              ierr= NF_INQ_DIMID(nid,"VI Wavenumber",id(3))
    241242              ierr= NF_INQ_DIMID(nid,"Time",id(4))
    242243
Note: See TracChangeset for help on using the changeset viewer.