Ignore:
Timestamp:
Mar 22, 2012, 8:25:03 PM (13 years ago)
Author:
jleconte
Message:
  • New turbulent diffusion scheme solving "most" energy conservation problems:
    • Turbulent energy created by buoyancy effects is now dissipated back into enthalpy
    • the scheme is now written in an enthalpy conservative way
  • Turbulent diffusion now treated in routine turbdiff (in F90).
  • Temporarily, for comparison, the old vdifc can be used

by setting UseTurbDiff?=.false. in physiq.F90

  • The sensible heat flux is now an output
  • Corrected evaporation at the surface when all the surface water is evaporated.
Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

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

    r588 r594  
    9898!           Improved water cycle: R. Wordsworth / B. Charnay (2010)
    9999!           To F90: R. Wordsworth (2010)
     100!           New turbulent diffusion scheme: J. Leconte (2012)
    100101!     
    101102!==================================================================
    … …  
    200201      real OLR_nu(ngridmx,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
    201202      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
     203 
     204      real sensibFlux(ngridmx)   ! turbulent flux given by the atm to the surface
     205 
     206      save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,sensibFlux
    203207
    204208
    … …  
    285289      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS
    286290
    287 !     included by RW for temporary comparison
    288       real zdtnirco2(ngridmx,nlayermx) ! (K/s)
    289 
    290291!     included by RW to compute tracer column densities
    291292      real qcol(ngridmx,nqmx)
    … …  
    311312      real dtsurfh2olat(ngridmx)
    312313
    313 !     included by RW to test energy conservation
    314       real dEtot, dEtots, masse, vabs, dvabs
     314!     to test energy conservation (RW+JL)
     315      real dEtot, dEtots, masse, AtmToSurf_TurbFlux
    315316      real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
    316 
     317      real dEzRad(ngridmx,nlayermx),dEzdif(ngridmx,nlayermx)
     318      real madjdE(ngridmx), lscaledE(ngridmx)
     319!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
     320     
    317321      real dItot, dVtot
    318 
    319 !     included by RW to test water conservation
    320       real h2otot
    321 
    322 
    323322
    324323!     included by BC for evaporation     
    … …  
    332331
    333332!     included by RW to test water conservation (by routine)
     333      real h2otot
    334334      real dWtot, dWtots
    335335      real h2o_surf_all
    … …  
    364364      real totcloudfrac(ngridmx)
    365365
    366 !     included by RW for Tmax test
    367       integer iTmax
    368 
    369366!     included by RW for vdifc water conservation test
    370367      real nconsMAX
    … …  
    392389      real reffcol(ngridmx,2)
    393390
    394 !     included by RW for flexible 3D energy conservation testing
    395       real vdifcdE(ngridmx), madjdE(ngridmx), lscaledE(ngridmx)
    396391
    397392!     included by RW for sourceevol
    … …  
    405400      logical ice_update
    406401      save ice_update
     402
     403!JL12 temporary to test difference in diffusion schemes
     404      logical UseTurbDiff
    407405
    408406!=======================================================================
    … …  
    555553
    556554            if(ngrid.eq.1)then
    557                qzero1D               = 0.0
     555               qzero1D               = 20.0
    558556               qsurf(1,igcm_h2o_vap) = qzero1D
    559557               do l=1, nlayer
    … …  
    642640         do ig=1,ngrid
    643641            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
    644             zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
    645642         enddo
    646643      enddo
    … …  
    839836            dEtotsLW = dEtotsLW/totarea
    840837            print*,'---------------------------------------------------------------'
    841             print*,'In corrk SW atmospheric energy change =',dEtotSW,' W m-2'
    842             print*,'In corrk LW atmospheric energy change =',dEtotLW,' W m-2'
    843             print*,'atmospheric energy change  (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
    844             print*,'In corrk SW surface energy change     =',dEtotsSW,' W m-2'
    845             print*,'In corrk LW surface energy change     =',dEtotsLW,' W m-2'
    846             print*,'surface energy change       (SW+LW)   =',dEtotsLW+dEtotsSW,' W m-2'
     838            print*,'In corrk SW atmospheric heating      =',dEtotSW,' W m-2'
     839            print*,'In corrk LW atmospheric heating      =',dEtotLW,' W m-2'
     840            print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
     841            print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
     842            print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
     843            print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
    847844         endif
    848845!-------------------------
    … …  
    862859         zdum1(:,:)=0.0
    863860         zdum2(:,:)=0.0
    864          do l=1,nlayer
    865             do ig=1,ngrid
     861
     862
     863!JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff     
     864         UseTurbDiff=.true.
     865         if (UseTurbDiff) then
     866         
     867           call turbdiff(ngrid,nlayer,nq,rnat,       &
     868              ptimestep,capcal,lwrite,               &
     869              pplay,pplev,zzlay,zzlev,z0,            &
     870              pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
     871              zdum1,zdum2,pdt,pdq,zflubid,           &
     872              zdudif,zdvdif,zdtdif,zdtsdif,          &
     873              sensibFlux,q2,zdqdif,zdqsdif,lastcall)
     874
     875         else
     876         
     877           do l=1,nlayer
     878             do ig=1,ngrid
     879               zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
    866880               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
    867             enddo
    868          enddo
    869 
    870          call vdifc(ngrid,nlayer,nq,rnat,zpopsk,     &
     881             enddo
     882           enddo
     883 
     884           call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
    871885              ptimestep,capcal,lwrite,               &
    872886              pplay,pplev,zzlay,zzlev,z0,            &
    873887              pu,pv,zh,pq,tsurf,emis,qsurf,          &
    874888              zdum1,zdum2,zdh,pdq,zflubid,           &
    875               zdudif,zdvdif,zdhdif,zdtsdif,q2,       &
    876               zdqdif,zdqsdif,lastcall)
     889              zdudif,zdvdif,zdhdif,zdtsdif,          &
     890              sensibFlux,q2,zdqdif,zdqsdif,lastcall)
     891
     892           do l=1,nlayer
     893             do ig=1,ngrid
     894                zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
     895             enddo
     896           enddo
     897
     898         end if !turbdiff
    877899
    878900         do l=1,nlayer
    … …  
    880902               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
    881903               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
    882                pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
    883                zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
     904               pdt(ig,l)=pdt(ig,l)+zdtdif(ig,l)
    884905            enddo
    885906         enddo
    … …  
    910931            dEtot=0.0
    911932            dEtots=0.0
    912 
    913             vdifcdE(:)=0.0
     933            dEzdif(:,:)=0.
     934            AtmToSurf_TurbFlux=0.
     935           
    914936            do ig = 1, ngrid
    915937               do l = 1, nlayer
    916938                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
    917939                  dEtot = dEtot + cpp*masse*zdtdif(ig,l)*area(ig)
    918                  
    919                   vabs  = sqrt(pdu(ig,l)**2 + pdv(ig,l)**2)
    920                   dvabs = sqrt(zdudif(ig,l)**2 + zdvdif(ig,l)**2)
    921                   dEtot = dEtot + masse*vabs*dvabs*area(ig)
    922 
    923                   vdifcdE(ig) = vdifcdE(ig) + masse*vabs*dvabs
    924 
    925                enddo
    926                dEtot = dEtot - zflubid(ig)*area(ig) ! subtract flux from ground
    927 
    928                dEtots = dEtots + capcal(ig)*zdtsdif(ig)*area(ig)
    929                vdifcdE(ig) = vdifcdE(ig) + capcal(ig)*zdtsdif(ig)
     940                  dEzdif (ig,l)= cpp*masse*(zdtdif(ig,l))*area(ig)
     941               enddo
     942               dEtot = dEtot + sensibFlux(ig)*area(ig)! subtract flux to the ground
     943               dEtots = dEtots + capcal(ig)*zdtsdif(ig)*area(ig)-zflubid(ig)*area(ig)-sensibFlux(ig)*area(ig)
     944               AtmToSurf_TurbFlux=AtmToSurf_TurbFlux+sensibFlux(ig)*area(ig)
    930945            enddo
    931946            dEtot  = dEtot/totarea
    932947            dEtots = dEtots/totarea
    933             print*,'In difv atmospheric energy change     =',dEtot,' W m-2'
    934             print*,'In difv surface energy change         =',dEtots,' W m-2'
    935             !print*,'Note we ignore the wind change...'
    936             ! and latent heat for that matter...
     948            AtmToSurf_TurbFlux=AtmToSurf_TurbFlux/totarea
     949            print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
     950            print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
     951            print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
     952! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
     953!    but not given back elsewhere
    937954         endif
    938955         !-------------------------
    … …  
    10531070            enddo
    10541071            dEtot=dEtot/totarea
    1055           print*,'In convadj atmospheric energy change    =',dEtot,' W m-2'
     1072          print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
    10561073         endif
    10571074         !-------------------------
    … …  
    16931710            OLR = OLR + area(ig)*fluxtop_lw(ig)
    16941711            GND = GND + area(ig)*fluxgrd(ig)
    1695             if(.not.callsoil) GND=GND+ area(ig)*fluxrad(ig)
    16961712            DYN = DYN + area(ig)*fluxdyn(ig)
    16971713           
    … …  
    17151731         
    17161732         if(enertest)then
    1717             print*,'SW energy balance SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
    1718             print*,'LW energy balance LW++ + ***ASR*** = ',dEtotLW+dEtotsLW+ASR/totarea,' W m-2'
    1719             print*,'LW energy balance LW++ - ***OLR*** = ',dEtotLW+dEtotsLW+OLR/totarea,' W m-2'
     1733            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
     1734            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR/totarea,' W m-2'
     1735            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR/totarea,' W m-2'
    17201736         endif
    17211737
    … …  
    19892005           call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
    19902006           call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
    1991            !call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
    1992            !call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
    1993            !call writediagfi(ngrid,"vdifcdE","heat from vdifc surface","W m-2",2,vdifcdE)
     2007        endif
     2008       
     2009        if(enertest) then
     2010          call writediagfi(ngridmx,"dEzdif","atm vdifc energy change","w.m^-2",3,dEzdif)
     2011          if(watercond) then
     2012             call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
     2013             call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
     2014          endif
    19942015        endif
    19952016
    … …  
    20012022
    20022023        ! debugging
    2003         !call writediagfi(ngrid,"vdifNC","H2O loss vdifc","kg m-2 s-1",2,vdifcncons)
    2004         !call writediagfi(ngrid,"cadjNC","H2O loss convadj","kg m-2 s-1",2,cadjncons)
    20052024        !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
    20062025        !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
    … …  
    20272046          endif
    20282047
    2029           if(watercond)then
     2048          if(watercond.or.CLFvarying)then
    20302049             !call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
    20312050             call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)

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

    r586 r594  
    44     &     pu,pv,ph,pq,ptsrf,pemis,pqsurf,
    55     &     pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
    6      &     pdudif,pdvdif,pdhdif,pdtsrf,pq2,
     6     &     pdudif,pdvdif,pdhdif,pdtsrf,sensibFlux,pq2,
    77     &     pdqdif,pdqsdif,lastcall)
    88
    … …  
    5555      REAL pfluxsrf(ngrid)
    5656      REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
    57       REAL pdtsrf(ngrid),pcapcal(ngrid)
     57      REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
    5858      REAL pq2(ngrid,nlay+1)
    5959     
    … …  
    674674         enddo
    675675      enddo
    676      
     676
     677      DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
     678         sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zh(ig,1)-ztsrf2(ig))
     679      ENDDO     
     680
    677681      if (tracer) then
    678682         do iq = 1, nq
Note: See TracChangeset for help on using the changeset viewer.