Changeset 3275 for trunk


Ignore:
Timestamp:
Mar 20, 2024, 3:05:14 PM (8 months ago)
Author:
afalco
Message:

Pluto PCM:
Changed _vap to _gas;
Included surfprop.F90;
callcorrk includes methane
AF

Location:
trunk
Files:
2 added
41 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.PLUTO.old

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/compile1d

    r3175 r3275  
    11#!/bin/bash
    2 
    3 ./makegcm_cygwin -debug -fdefault-real-8 -d 25 -b 17x23 -t 7 -s 1 -p pluto testphys1d
     2export PATH=$PATH:.
     3./makegcm_spirit_gfortran -debug -fdefault-real-8 -d 25 -t 7 -s 2 -b 17x23 -p pluto testphys1d
     4# ./makegcm_cygwin -debug -fdefault-real-8 -d 25 -b 17x23 -t 7 -s 1 -p pluto testphys1d

  • trunk/LMDZ.PLUTO.old/deftank

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/deftank/gcm

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/deftank/kbo_def

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/deftank/nogcm

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/deftank/nogcm_simple

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/deftank/testphys1d

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/bibio

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/dyn3d

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/dyn3d/poubelle

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/dyn3d/startsHD

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/dyn3d/stock

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/filtrez

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/grid

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/grid/dimension

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/phypluto

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/phypluto/callcorrk.F

    r3175 r3275  
    55     &     fluxsurf_sw,fluxtop_lw,fluxtop_sw,fluxtop_dn,
    66     &     reffrad,tau_col,ptime,pday,firstcall,lastcall,zzlay)
    7          
     7
    88      use radinc_h
    99      use radcommon_h
    10       use ioipsl_getincom 
     10      use ioipsl_getincom
    1111      use radii_mod
    1212      use aerosol_mod
    … …  
    2424!
    2525!     Authors
    26 !     ------- 
     26!     -------
    2727!     Emmanuel 01/2001, Forget 09/2001
    2828!     Robin Wordsworth (2009)
    … …  
    3737!-----------------------------------------------------------------------
    3838!     Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
    39 !     Layer #1 is the layer near the ground. 
     39!     Layer #1 is the layer near the ground.
    4040!     Layer #nlayermx is the layer at the top.
    4141
    … …  
    4444      INTEGER ngrid,nlayer
    4545      INTEGER igout
    46       REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol opacity tau 
     46      REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol opacity tau
    4747      REAL albedo(ngrid)                    ! SW albedo
    4848      REAL emis(ngrid)                      ! LW emissivity
    … …  
    110110      real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
    111111      REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
    112       real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) ! 
    113       real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) ! 
     112      real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) !
     113      real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) !
    114114      REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
    115115      REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
    … …  
    137137      save qxvaer, qsvaer, gvaer
    138138      save qxiaer, qsiaer, giaer
    139       save QREFvis3d, QREFir3d 
     139      save QREFvis3d, QREFir3d
    140140
    141141      REAL tau_col(ngrid) ! diagnostic from aeropacity
    … …  
    220220         call setspv            ! basic visible properties
    221221
    222          ! Radiative Hazes 
     222         ! Radiative Hazes
    223223         if (aerohaze) then
    224224
    … …  
    231231           !--------------------------------------------------
    232232           do iaer=1,naerkind
    233               if ((iaer.eq.iaero_haze)) then 
    234                call haze_reffrad(ngrid,nlayer,reffrad(1,1,iaer), 
     233              if ((iaer.eq.iaero_haze)) then
     234               call haze_reffrad(ngrid,nlayer,reffrad(1,1,iaer),
    235235     &             nueffrad(1,1,iaer))
    236236              endif
    … …  
    238238           if (haze_radproffix) then
    239239              print*, 'haze_radproffix=T : fixed profile for haze rad'
    240            else 
     240           else
    241241              print*,'reffrad haze:',reffrad(1,1,iaero_haze)
    242242              print*,'nueff haze',nueffrad(1,1,iaero_haze)
    … …  
    272272!-----------------------------------------------------------------------
    273273!     Get 3D aerosol optical properties.
    274       ! ici on selectionne les proprietes opt correspondant a reffrad 
     274      ! ici on selectionne les proprietes opt correspondant a reffrad
    275275      if (aerohaze) then
    276276        !--------------------------------------------------
    … …  
    283283        endif
    284284
    285         call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,         
    286      &       QVISsQREF3d,omegaVIS3d,gVIS3d,                         
    287      &      QIRsQREF3d,omegaIR3d,gIR3d,                             
     285        call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
     286     &       QVISsQREF3d,omegaVIS3d,gVIS3d,
     287     &      QIRsQREF3d,omegaIR3d,gIR3d,
    288288     &      QREFvis3d,QREFir3d)
    289289
    290290        ! Get aerosol optical depths.
    291         call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,     
    292      &      reffrad,QREFvis3d,QREFir3d,                             
     291        call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,
     292     &      reffrad,QREFvis3d,QREFir3d,
    293293     &      tau_col)
    294294      endif
    … …  
    298298      IF (methane) then
    299299        vmrch4(:,:)=0.
    300        
     300
    301301        if (ch4fix) then
    302302           if (vmrch4_proffix) then
    303303            !! Interpolate on the model vertical grid
    304304             do ig=1,ngridmx
    305                CALL interp_line(levdat,vmrdat,Nfine,
    306     &                  zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
     305    !            CALL interp_line(levdat,vmrdat,Nfine,
     306    ! &                  zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
    307307             enddo
    308308           else
    … …  
    317317!     Prepare NON LTE correction in Pluto atmosphere
    318318      IF (nlte) then
    319         CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,
    320     &             eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
     319    !     CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,
     320    ! &             eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
    321321      ENDIF
    322322c     Net atmospheric radiative cooling rate from C2H2 (K.s-1):
    … …  
    343343!     shortwave
    344344            do iaer=1,naerkind
    345                DO nw=1,L_NSPECTV 
     345               DO nw=1,L_NSPECTV
    346346                  do l=1,nlayermx
    347347
    348                      temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer) 
     348                     temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer)
    349349     $                    *QREFvis3d(ig,nlayermx+1-l,iaer)
    350350
    351                      temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer) 
     351                     temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer)
    352352     $                    *QREFvis3d(ig,max(nlayermx-l,1),iaer)
    353353                     qxvaer(2*l,nw,iaer)  = temp1
    … …  
    378378
    379379!     longwave
    380                DO nw=1,L_NSPECTI 
     380               DO nw=1,L_NSPECTI
    381381                  do l=1,nlayermx
    382382
    383                      temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer) 
     383                     temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer)
    384384     $                    *QREFir3d(ig,nlayermx+1-l,iaer)
    385385
    386                      temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer) 
     386                     temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer)
    387387     $                    *QREFir3d(ig,max(nlayermx-l,1),iaer)
    388388
    … …  
    421421                  do nw=1,L_NSPECTV
    422422                     if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
    423                         print*,'Serious problems with qsvaer values' 
     423                        print*,'Serious problems with qsvaer values'
    424424                        print*,'in callcorrk'
    425425                        call abort
    … …  
    430430                  end do
    431431
    432                   do nw=1,L_NSPECTI 
     432                  do nw=1,L_NSPECTI
    433433                     if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
    434434                        print*,'Serious problems with qsiaer values'
    … …  
    448448!-----------------------------------------------------------------------
    449449!     Aerosol optical depths
    450          IF (aerohaze) THEN   
    451           do iaer=1,naerkind     ! heritage generic       
     450         IF (aerohaze) THEN
     451          do iaer=1,naerkind     ! heritage generic
    452452            do k=0,nlayer-1
    453453               pweight=
    454454     $              (pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/
    455455     $              (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
    456                if (QREFvis3d(ig,L_NLAYRAD-k,iaer).ne.0) then 
     456               if (QREFvis3d(ig,L_NLAYRAD-k,iaer).ne.0) then
    457457                 temp=aerosol(ig,L_NLAYRAD-k,iaer)/
    458458     $              QREFvis3d(ig,L_NLAYRAD-k,iaer)
    … …  
    479479!     Albedo and emissivity
    480480         albi=1-emis(ig)        ! longwave
    481          albv=albedo(ig)        ! shortwave 
     481         albv=albedo(ig)        ! shortwave
    482482         acosz=mu0(ig)          ! cosine of sun incident angle
    483483
    484484!-----------------------------------------------------------------------
    485 !     Methane vapour 
     485!     Methane vapour
    486486
    487487c     qvar = mixing ratio
    … …  
    490490c     datagcm/composition.in for the k-coefficients.
    491491         qvar(:)=0.
    492          IF (methane) then 
     492         IF (methane) then
    493493
    494494           do l=1,nlayer
    … …  
    554554!! following lines changed in 03/2015 to solve upper atmosphere bug
    555555!        plevrad(1) = 0.
    556 !        plevrad(2) = max(pgasmin,0.0001*plevrad(3))     
     556!        plevrad(2) = max(pgasmin,0.0001*plevrad(3))
    557557!
    558558!        tlevrad(1) = tlevrad(2)
    … …  
    563563!
    564564!        pmid(1) = plevrad(2)
    565 !        pmid(2) = plevrad(2) 
     565!        pmid(2) = plevrad(2)
    566566
    567567         DO l=1,L_NLAYRAD-1
    … …  
    574574         pmid(L_LEVELS) = plevrad(L_LEVELS)
    575575         tmid(L_LEVELS) = tlevrad(L_LEVELS)
    576      
     576
    577577      !TB
    578578         if ((PMID(2).le.1.e-5).and.(ig.eq.1)) then
    … …  
    608608          endif
    609609         enddo
    610      
     610
    611611!=======================================================================
    612612!     Calling the main radiative transfer subroutines
    … …  
    614614!-----------------------------------------------------------------------
    615615!     Shortwave
    616        
     616
    617617         IF(fract(ig) .GE. 1.0e-4) THEN ! only during daylight  IPM?! flux UV...
    618618
    … …  
    623623            END DO
    624624
    625             !print*, 'starting optcv' 
     625            !print*, 'starting optcv'
    626626            call optcv(dtauv,tauv,taucumv,plevrad,
    627627     $           qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,
    … …  
    666666!        IR spectral output from top of the atmosphere
    667667         if(specOLR)then
    668             do nw=1,L_NSPECTI 
     668            do nw=1,L_NSPECTI
    669669               OLR_nu(ig,nw)=nfluxtopi_nu(nw)
    670670            end do
    … …  
    673673! **********************************************************
    674674!     Finally, the heating rates
    675 !     g/cp*DF/DP 
     675!     g/cp*DF/DP
    676676! **********************************************************
    677677
    … …  
    682682               !dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))*dpp   !averaged dtlw on each wavelength
    683683               do nw=1,L_NSPECTV
    684                  dtsw_nu(L_NLAYRAD+1-l,nw)= 
     684                 dtsw_nu(L_NLAYRAD+1-l,nw)=
    685685     &              (fmnetv_nu(l,nw)-fmnetv_nu(l-1,nw))*dpp
    686686               end do
    … …  
    689689               !dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))*dpp   !averaged dtlw on each wavelength
    690690               do nw=1,L_NSPECTI
    691                  dtlw_nu(L_NLAYRAD+1-l,nw)= 
    692      &              (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp     
     691                 dtlw_nu(L_NLAYRAD+1-l,nw)=
     692     &              (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp
    693693               end do
    694          END DO     
    695          
     694         END DO
     695
    696696         ! values at top of atmosphere
    697697         dpp = g/(cpp*scalep*(plevrad(3)-plevrad(1)))
    698698
    699          ! SW 
     699         ! SW
    700700         !dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)*dpp
    701701         do nw=1,L_NSPECTV
    … …  
    704704         end do
    705705
    706          ! LW 
    707 c        dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) *dpp 
     706         ! LW
     707c        dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) *dpp
    708708         do nw=1,L_NSPECTI
    709709             dtlw_nu(L_NLAYRAD,nw)=
    … …  
    717717
    718718         if (.not.nlte) then
    719             eps_nlte_sw23(ig,:) =1. ! IF no NLTE 
    720             eps_nlte_sw33(ig,:) =1. ! IF no NLTE 
    721             eps_nlte_lw(ig,:) =1. ! IF no NLTE 
     719            eps_nlte_sw23(ig,:) =1. ! IF no NLTE
     720            eps_nlte_sw33(ig,:) =1. ! IF no NLTE
     721            eps_nlte_lw(ig,:) =1. ! IF no NLTE
    722722         endif
    723  
     723
    724724         do l=1,nlayer
    725725
    726726            !LW
    727             dtlw(ig,l) =0. 
     727            dtlw(ig,l) =0.
    728728!           dtlw_co(ig,l) =0.  ! only for diagnostic
    729729            do nw=1,L_NSPECTI
    730730              ! wewei : wavelength in micrometer
    731               if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then 
     731              if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then
    732732                dtlw_nu(l,nw)=dtlw_nu(l,nw)*eps_nlte_lw(ig,l)
    733               else 
     733              else
    734734                !dtlw_nu(l,nw)=1.*dtlw_nu(l,nw) ! no CO correction (Strobbel 1996)
    735735                dtlw_nu(l,nw)=0.33*dtlw_nu(l,nw) ! CO correction (Strobbel 1996)
    736736!               dtlw_co(ig,l)=dtlw_co(ig,l)+ dtlw_nu(l,nw) ! diagnostic
    737737              end if
    738               dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength 
     738              dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength
    739739            end do
    740740            ! adding c2h2 if cooling active
    … …  
    743743            !SW
    744744            dtsw(ig,l) =0.
    745  
     745
    746746            if (strobel) then
    747747
    748748             do nw=1,L_NSPECTV
    749               if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then 
     749              if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then
    750750                dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
    751               elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then 
     751              elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then
    752752                dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw33(ig,l)
    753753              else
    754754                dtsw_nu(l,nw)=dtsw_nu(l,nw)
    755755              end if
    756               dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw) 
     756              dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
    757757             end do
    758758
    … …  
    764764             enddo
    765765
    766             endif 
     766            endif
    767767
    768768
    … …  
    771771
    772772!     Diagnotics for last call for each grid point
    773          !if (lastcall) then 
     773         !if (lastcall) then
    774774
    775775          !print*,'albedi vis=',albv
    … …  
    806806      endif
    807807
    808       if(lastcall)then 
     808      if(lastcall)then
    809809
    810810        ! 1D Output
    … …  
    816816               open(116,file='surf_vals.out')
    817817               write(116,*) tsurf(1),pplev(1,1),
    818      &             fluxtop_dn(1) - fluxtop_sw(1),fluxtop_lw(1) 
     818     &             fluxtop_dn(1) - fluxtop_sw(1),fluxtop_lw(1)
    819819               do nw=1,L_NSPECTV
    820820                 write(116,*) wavev(nw),fmnetv_nu(L_NLAYRAD,nw)
    … …  
    830830           if(diagrad_OLR)then
    831831               open(117,file='OLRnu.out')
    832                write(117,*) 'IR wavel - band width - OLR' 
     832               write(117,*) 'IR wavel - band width - OLR'
    833833               do nw=1,L_NSPECTI
    834834                   write(117,*) wavei(nw),
    835835     &        abs(1.e4/bwnv(nw)-1.e4/bwnv(nw+1)),OLR_nu(1,nw)
    836                enddo   
     836               enddo
    837837               close(117)
    838838           endif
    … …  
    846846                 write(118,*) plevrad(2*l)
    847847                 do nw=1,L_NSPECTI
    848                      write(119,*) fluxupi_nu(l,nw) 
     848                     write(119,*) fluxupi_nu(l,nw)
    849849                 enddo
    850                enddo 
     850               enddo
    851851               close(118)
    852852               close(119)

  • trunk/LMDZ.PLUTO.old/libf/phypluto/cooling_stock

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO.old/libf/phypluto/physiq.F

    r3237 r3275  
    18491849          call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
    18501850          call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
    1851          ! call WRITEDIAGFI(ngrid,"pressure","Pression","Pa",3,pplay)
     1851          call WRITEDIAGFI(ngrid,"p","Pression","Pa",3,pplay)
    18521852          call WRITEDIAGFI(ngrid,"fluxrad","fluxrad",
    18531853     &                                        "W m-2",2,fluxrad)
    … …  
    20732073          call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",0,ps)
    20742074          call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
     2075          call WRITEDIAGFI(ngrid,"p","Pression","Pa",3,pplay)
    20752076
    20762077          call WRITEDIAGFI(ngrid,"fluxsurf_sw","sw surface flux",

  • trunk/LMDZ.PLUTO.old/libo

    • Property svn:ignore set to
      *.e
      *.mod

  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/lect_start_archive.F

    r3198 r3275  
    776776      DO iq=1,nqtot
    777777          txt=trim(tname(iq))//"_surf"
    778           if (txt.eq."h2o_vap") then
    779             ! There is no surface tracer for h2o_vap;
     778          if (txt.eq."h2o_gas") then
     779            ! There is no surface tracer for h2o_gas;
    780780            ! "h2o_ice" should be loaded instead
    781781            txt="h2o_ice_surf"
    782782            write(*,*) 'lect_start_archive: loading surface tracer',
    783      &                     ' h2o_ice instead of h2o_vap'
     783     &                     ' h2o_ice instead of h2o_gas'
    784784          endif
    785785

  • trunk/LMDZ.PLUTO/libf/phypluto

    • Property svn:ignore set to
      *.mod
      vdifc_mod.F.save

  • trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90

    r3195 r3275  
    9494
    9595      real CLFtot
    96       integer igen_ice,igen_vap ! to store the index of generic tracer
     96      integer igen_ice,igen_gas ! to store the index of generic tracer
    9797      logical dummy_bool ! dummy boolean just in case we need one
    9898      ! integer i_rgcs_ice(aerogeneric)

  • trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90

    r3197 r3275  
    77      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
    88          albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    &
    9           tsurf,fract,dist_star,aerosol,muvar,                 &
     9          zzlay,tsurf,fract,dist_star,aerosol,muvar,                 &
    1010          dtlw,dtsw,fluxsurf_lw,                               &
    1111          fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
    … …  
    3030      use aeropacity_mod, only: aeropacity
    3131      use aeroptproperties_mod, only: aeroptproperties
    32       use tracer_h, only: constants_epsi_generic
    33       use comcstfi_mod, only: pi, mugaz, cpp
     32      use tracer_h, only: constants_epsi_generic,igcm_ch4_gas,igcm_n2,mmol
     33      use comcstfi_mod, only: pi, mugaz, cpp, r, g
    3434      use callkeys_mod, only: varactive,diurnal,tracer,varfixed,satval, &
    3535                              diagdtau,kastprof,strictboundcorrk,specOLR, &
    3636                              tplanckmin,tplanckmax,global1d, &
    37                               generic_condensation, aerohaze, haze_radproffix
     37                              generic_condensation,aerohaze,haze_radproffix,&
     38                              methane,carbox,cooling,nlte,strobel,&
     39                              ch4fix,vmrch4_proffix,vmrch4fix
    3840      use optcv_mod, only: optcv
    3941      use optci_mod, only: optci
    … …  
    4446      use generic_tracer_index_mod, only: generic_tracer_index
    4547      use planetwide_mod, only: planetwide_maxval, planetwide_minval
     48      use radcommon_h, only: wavev,wavei
    4649      implicit none
    4750
    … …  
    8083      REAL,INTENT(IN) :: pplay(ngrid,nlayer)       ! Mid-layer pressure (Pa).
    8184      REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
     85      REAL,INTENT(IN) :: zzlay(ngrid,nlayer)       ! Mid-layer altitude
    8286      REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
    8387      REAL,INTENT(IN) :: fract(ngrid)              ! Fraction of day.
    … …  
    157161!$OMP THREADPRIVATE(tauaero)
    158162      REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
     163      real*8 nfluxtopv_nu(L_NSPECTV)
    159164      REAL*8 nfluxoutv_nu(L_NSPECTV)                 ! Outgoing band-resolved VI flux at TOA (W/m2).
    160165      REAL*8 nfluxtopi_nu(L_NSPECTI)                 ! Net band-resolved IR flux at TOA (W/m2).
    161166      REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)         ! For 1D diagnostic.
    162167      REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
     168      real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) !
     169      real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) !
    163170      REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
    164171      REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
    … …  
    211218      REAL*8 surface_stellar_flux   ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
    212219
     220      !     NLTE factor for CH4
     221      real eps_nlte_sw23(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw
     222      real eps_nlte_sw33(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw
     223      real eps_nlte_lw(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw
     224      integer Nfine,ifine
     225      parameter(Nfine=701)
     226      real,save :: levdat(Nfine),vmrdat(Nfine)
     227      real :: vmrch4(ngrid,nlayer)              ! vmr ch4 from vmrch4_proffix
     228
     229      REAL dtlw_nu(nlayer,L_NSPECTI) ! heating rate (K/s) due to LW in spectral bands
     230      REAL dtsw_nu(nlayer,L_NSPECTV) ! heating rate (K/s) due to SW in spectral bands
     231
    213232      ! local variable
     233      REAL dpp  ! intermediate
     234
    214235      integer ok ! status (returned by NetCDF functions)
    215236
    216       integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
    217       logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
     237      integer igcm_generic_gas, igcm_generic_ice! index of the vap and ice of generic_tracer
     238      logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer
    218239      real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
    219240!$OMP THREADPRIVATE(metallicity)
    … …  
    327348         ENDIF
    328349
    329 #ifndef MESOSCALE
    330350         if (is_master) call system('rm -f surf_vals_long.out')
    331 #endif
    332351
    333352         call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
    … …  
    434453
    435454
    436          if(varfixed .and. generic_condensation)then
     455         if(varfixed .and. (generic_condensation .or. methane .or. carbox))then
    437456            write(*,*) "Deep generic tracer vapor mixing ratio ? (no effect if negative) "
    438457            qvap_deep=-1. ! default value
    … …  
    504523      endif
    505524
     525
     526!-----------------------------------------------------------------------
     527!     Prepare CH4 mixing ratio for radiative transfer
     528      IF (methane) then
     529        vmrch4(:,:)=0.
     530
     531        if (ch4fix) then
     532           if (vmrch4_proffix) then
     533            !! Interpolate on the model vertical grid
     534             do ig=1,ngrid
     535               CALL interp_line(levdat,vmrdat,Nfine, &
     536                      zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
     537             enddo
     538           else
     539            vmrch4(:,:)=vmrch4fix
     540           endif
     541        else
     542            vmrch4(:,:)=pq(:,:,igcm_ch4_gas)*100.* &
     543                        mmol(igcm_n2)/mmol(igcm_ch4_gas)
     544        endif
     545      ENDIF
     546
     547!     Prepare NON LTE correction in Pluto atmosphere
     548      IF (nlte) then
     549        CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,&
     550                  eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
     551      ENDIF
     552!     Net atmospheric radiative cooling rate from C2H2 (K.s-1):
     553!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     554    !   dtlw_hcn_c2h2=0.
     555      if (cooling) then
     556        !  CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,&
     557                                !   pt,dtlw_hcn_c2h2)
     558      endif
     559
     560
    506561!-----------------------------------------------------------------------
    507562      do ig=1,ngrid ! Starting Big Loop over every GCM column
    … …  
    671726      !-----------------------------------------------------------------------
    672727
    673       if (generic_condensation) then
     728      if (generic_condensation .or. methane .or. carbox) then
    674729
    675730         ! For now, only one GCS tracer can be both variable and radiatively active
    … …  
    679734         do iq=1,nq
    680735
    681             call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    682 
    683             if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
     736            call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
     737
     738            if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
    684739
    685740               if(varactive)then
    686741
    687                   i_var=igcm_generic_vap
     742                  i_var=igcm_generic_gas
    688743                  do l=1,nlayer
    689744                     qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
    … …  
    731786      !-----------------------------------------------------------------------
    732787
    733       if (.not. generic_condensation) then
     788      if (.not. (generic_condensation .or. methane .or. carbox)) then
    734789         do k=1,L_LEVELS
    735790            qvar(k) = 1.0D-7
    … …  
    741796         ! IMPORTANT: Now convert from kg/kg to mol/mol.
    742797         do k=1,L_LEVELS
    743             if (generic_condensation) then
     798            if (generic_condensation .or. methane .or. carbox) then
    744799               do iq=1,nq
    745                   call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    746                   if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
    747 
    748                      epsi_generic=constants_epsi_generic(iq)
    749 
    750                      qvar(k) = qvar(k)/(epsi_generic+qvar(k)*(1.-epsi_generic))
     800                  call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
     801                  if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
     802                    if(.not. varactive .or. i_var.eq.iq)then
     803
     804                        epsi_generic=constants_epsi_generic(iq)
     805
     806                        qvar(k) = qvar(k)/(epsi_generic+qvar(k)*(1.-epsi_generic))
     807                    endif
    751808
    752809                  endif
    … …  
    762819      DO l=1,nlayer
    763820         muvarrad(2*l)   = muvar(ig,nlayer+2-l)
    764          muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
     821             muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
    765822      END DO
    766823
    … …  
    9651022                 tmid,pmid,taugsurf,qvar,muvarrad)
    9661023
    967             call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
    968                  acosz,stel_fract,                                 &
    969                  nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu,  &
    970                  fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
     1024            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,   &
     1025                 acosz,stel_fract,nfluxtopv,fluxtopvdn,nfluxoutv_nu,&
     1026                 nfluxgndv_nu,nfluxtopv_nu, &
     1027                 fmnetv,fmnetv_nu,fluxupv,fluxdnv,fzerov,taugsurf)
    9711028
    9721029         else ! During the night, fluxes = 0.
    … …  
    10101067         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
    10111068              wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu,         &
    1012               fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
     1069              fmneti,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
    10131070
    10141071!-----------------------------------------------------------------------
    … …  
    10501107
    10511108!     Finally, the heating rates
    1052 
    10531109         DO l=2,L_NLAYRAD
    1054             dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
    1055                 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
    1056             dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
    1057                 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
     1110            ! dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
     1111            !     *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
     1112            dpp = glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
     1113            do nw=1,L_NSPECTV
     1114                dtsw_nu(L_NLAYRAD+1-l,nw)= &
     1115                   (fmnetv_nu(l,nw)-fmnetv_nu(l-1,nw))*dpp
     1116            end do
     1117
     1118            ! dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
     1119            !     *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
     1120            do nw=1,L_NSPECTI
     1121                dtlw_nu(L_NLAYRAD+1-l,nw)= &
     1122                   (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp
     1123            end do
    10581124         END DO
    10591125
    10601126!     These are values at top of atmosphere
    1061          dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
    1062              *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
    1063          dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
    1064              *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
     1127        !  dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
     1128        !      *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
     1129        !  dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
     1130        !      *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
     1131         dpp = g/(cpp*scalep*(plevrad(3)-plevrad(1)))
     1132         do nw=1,L_NSPECTV
     1133            dtsw_nu(L_NLAYRAD,nw)= &
     1134             (fmnetv_nu(1,nw)-nfluxtopv_nu(nw))*dpp
     1135         end do
     1136         do nw=1,L_NSPECTI
     1137             dtlw_nu(L_NLAYRAD,nw)= &
     1138             (fmneti_nu(1,nw)-nfluxtopi_nu(nw))*dpp
     1139         end do
    10651140
    10661141      !  Optical thickness diagnostics (added by JVO)
    … …  
    10841159      endif
    10851160
     1161! **********************************************************
     1162!     NON NLTE correction in Pluto atmosphere
     1163!     And conversion of LW spectral heating rates to total rates
     1164! **********************************************************
     1165
     1166      if (.not.nlte) then
     1167        eps_nlte_sw23(ig,:) =1. ! IF no NLTE
     1168        eps_nlte_sw33(ig,:) =1. ! IF no NLTE
     1169        eps_nlte_lw(ig,:) =1. ! IF no NLTE
     1170     endif
     1171
     1172     do l=1,nlayer
     1173
     1174        !LW
     1175        dtlw(ig,l) =0.
     1176!           dtlw_co(ig,l) =0.  ! only for diagnostic
     1177        do nw=1,L_NSPECTI
     1178          ! wewei : wavelength in micrometer
     1179          if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then
     1180            dtlw_nu(l,nw)=dtlw_nu(l,nw)*eps_nlte_lw(ig,l)
     1181          else
     1182            !dtlw_nu(l,nw)=1.*dtlw_nu(l,nw) ! no CO correction (Strobbel 1996)
     1183            dtlw_nu(l,nw)=0.33*dtlw_nu(l,nw) ! CO correction (Strobbel 1996)
     1184!               dtlw_co(ig,l)=dtlw_co(ig,l)+ dtlw_nu(l,nw) ! diagnostic
     1185          end if
     1186          dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength
     1187        end do
     1188        ! adding c2h2 if cooling active
     1189        ! dtlw(ig,l)=dtlw(ig,l)+dtlw_hcn_c2h2(ig,l)
     1190
     1191        !SW
     1192        dtsw(ig,l) =0.
     1193
     1194        if (strobel) then
     1195
     1196         do nw=1,L_NSPECTV
     1197          if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then
     1198            dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
     1199          elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then
     1200            dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw33(ig,l)
     1201          else
     1202            dtsw_nu(l,nw)=dtsw_nu(l,nw)
     1203          end if
     1204          dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
     1205         end do
     1206
     1207        else ! total heating rates multiplied by nlte coef
     1208
     1209         do nw=1,L_NSPECTV
     1210            dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
     1211            dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
     1212         enddo
     1213
     1214        endif
     1215
     1216
     1217     end do
     1218! **********************************************************
     1219
    10861220
    10871221!-----------------------------------------------------------------------
    10881222      end do ! End of big loop over every GCM column.
    10891223!-----------------------------------------------------------------------
    1090 
    10911224
    10921225

  • trunk/LMDZ.PLUTO/libf/phypluto/condensation_generic_mod.F90

    r3184 r3275  
    1818!     -------
    1919!     Calculates large-scale condensation of generic tracer "tname".
    20 !     By convention, tname ends with the suffix "_vap", as it represents the
     20!     By convention, tname ends with the suffix "_gas", as it represents the
    2121!     gas phase of the generic tracer
    2222!     
    … …  
    3939        REAL, intent(in) :: pdt(ngrid,nlayer)     ! physical temperature tendency (K/s)
    4040        REAL, intent(in) :: pdq(ngrid,nlayer,nq)  ! physical tracer tendency (K/s)
    41         ! CHARACTER(*), intent(in) :: tname_vap     ! name of the tracer we consider. BY convention, it ends with _vap !!!
     41        ! CHARACTER(*), intent(in) :: tname_gas     ! name of the tracer we consider. BY convention, it ends with _gas !!!
    4242        REAL, intent(out) :: pdtlsc(ngrid,nlayer)  ! incrementation de la temperature (K)
    4343        REAL, intent(out) :: pdqvaplsc(ngrid,nlayer,nq) ! incrementation de la vapeur du traceur
    … …  
    4949        REAL, SAVE :: qvap_deep   ! deep mixing ratio of vapor when simulating bottom less planets
    5050        REAL, SAVE :: qvap_top   ! top mixing ratio of vapor when simulating bottom less planets
    51         logical, save :: perfect_vap_profile
    52 !$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, perfect_vap_profile)
     51        logical, save :: perfect_gas_profile
     52!$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, perfect_gas_profile)
    5353
    5454!       Local variables
    5555
    5656        ! to call only one time the ice/vap pair of a tracer
    57         logical call_ice_vap_generic
     57        logical call_ice_gas_generic
    5858
    5959        INTEGER i, k , nn, iq
    … …  
    6969        real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,Lcp,zqs_temp,zdqs
    7070        real zqs_temp_1, zqs_temp_2, zqs_temp_3
    71         integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
     71        integer igcm_generic_gas, igcm_generic_ice! index of the vap and ice of generic_tracer
    7272        ! CHARACTER(len=*) :: tname_ice
    7373        ! evaporation calculations
    … …  
    9696                write(*,*) " qvap_top = ",qvap_top
    9797               
    98                 write(*,*) " perfect_vap_profile ? "
    99                 perfect_vap_profile=.false. ! default value
    100                 call getin_p("perfect_vap_profile",perfect_vap_profile)
    101                 write(*,*) " perfect_vap_profile = ",perfect_vap_profile
     98                write(*,*) " perfect_gas_profile ? "
     99                perfect_gas_profile=.false. ! default value
     100                call getin_p("perfect_gas_profile",perfect_gas_profile)
     101                write(*,*) " perfect_gas_profile = ",perfect_gas_profile
    102102
    103103                firstcall = .false.
    … …  
    115115        do iq=1,nq
    116116
    117                 call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    118 
    119                 if(call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
     117                call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
     118
     119                if(call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
    120120                        m=constants_mass(iq)
    121                         delta_vapH=constants_delta_vapH(iq)
     121                        delta_gasH=constants_delta_gasH(iq)
    122122                        Tref=constants_Tref(iq)
    123123                        Pref=constants_Pref(iq)
    … …  
    139139                                        !           zt(i)=15.   ! check too low temperatures
    140140                                        endif
    141                                         zx_q(i) = pq(i,k,igcm_generic_vap)+pdq(i,k,igcm_generic_vap)*ptimestep
     141                                        zx_q(i) = pq(i,k,igcm_generic_gas)+pdq(i,k,igcm_generic_gas)*ptimestep
    142142
    143143                                        call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp)
    … …  
    189189
    190190                                !Tendances de t et q
    191                                 pdqvaplsc(1:ngrid,k,igcm_generic_vap)  = - zcond(1:ngrid)
    192                                 pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_vap)
     191                                pdqvaplsc(1:ngrid,k,igcm_generic_gas)  = - zcond(1:ngrid)
     192                                pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_gas)
    193193                                pdtlsc(1:ngrid,k)  = pdtlsc(1:ngrid,k) + pdqliqlsc(1:ngrid,k,igcm_generic_ice)*Lcp
    194194
    195195                        Enddo ! k= nlayer, 1, -1
    196196
    197                         if ((perfect_vap_profile) .and. (ngrid.eq.1)) then
    198                                 ! perfect_vap_profile is a mode that should a priori only be used in 1D:
     197                        if ((perfect_gas_profile) .and. (ngrid.eq.1)) then
     198                                ! perfect_gas_profile is a mode that should a priori only be used in 1D:
    199199                                ! as it is written below, it aims to force the vap profile:
    200200                                ! - below condensation, it is forced to qvap_deep
    … …  
    202202                                ! - above the cold trap, it is forced to the value allowed by the cold trap
    203203
    204                                 ! perfect_vap_profile can be customed as you want
     204                                ! perfect_gas_profile can be customed as you want
    205205                               
    206206                                tau = 10. * ptimestep ! tau must not be lower than the physical step time.
    … …  
    223223                                if (k_cold_trap .lt. nlayer) then
    224224                                        do k=k_cold_trap+1,nlayer
    225                                                 pdqvaplsc(1,k,igcm_generic_vap) = (pq(1,k_cold_trap,igcm_generic_vap) - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)
     225                                                pdqvaplsc(1,k,igcm_generic_gas) = (pq(1,k_cold_trap,igcm_generic_gas) - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas)
    226226                                        enddo
    227227                                endif
    … …  
    231231                                        call Psat_generic(zt(1),pplay(1,k),metallicity,psat_tmp,zqs_temp)
    232232                                        if (zqs_temp .gt. qvap_deep) then
    233                                                 pdqvaplsc(1,k,igcm_generic_vap)  = (qvap_deep - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)
     233                                                pdqvaplsc(1,k,igcm_generic_gas)  = (qvap_deep - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas)
    234234                                        endif
    235235                                        if (zqs_temp .lt. qvap_deep) then
    236                                                 pdqvaplsc(1,k,igcm_generic_vap)  = (0.99*zqs_temp - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)
     236                                                pdqvaplsc(1,k,igcm_generic_gas)  = (0.99*zqs_temp - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas)
    237237                                        endif
    238238                                enddo
    … …  
    246246                                ! brings lower generic vapor ratio to a fixed value.
    247247                                ! tau is in seconds and must not be lower than the physical step time.
    248                                 pdqvaplsc(1:ngrid,1,igcm_generic_vap) = (qvap_deep - pq(1:ngrid,1,igcm_generic_vap))/tau - pdq(1:ngrid,1,igcm_generic_vap)
     248                                pdqvaplsc(1:ngrid,1,igcm_generic_gas) = (qvap_deep - pq(1:ngrid,1,igcm_generic_gas))/tau - pdq(1:ngrid,1,igcm_generic_gas)
    249249                        endif
    250250                        if (qvap_top >= 0.) then
    … …  
    252252                                ! brings lower generic vapor ratio to a fixed value.
    253253                                ! tau is in seconds and must not be lower than the physical step time.
    254                                 pdqvaplsc(1:ngrid,nlayer,igcm_generic_vap) = (qvap_top - pq(1:ngrid,nlayer,igcm_generic_vap))/tau - pdq(1:ngrid,nlayer,igcm_generic_vap)
     254                                pdqvaplsc(1:ngrid,nlayer,igcm_generic_gas) = (qvap_top - pq(1:ngrid,nlayer,igcm_generic_gas))/tau - pdq(1:ngrid,nlayer,igcm_generic_gas)
    255255                        endif
    256                 endif !if(call_ice_vap_generic)
     256                endif !if(call_ice_gas_generic)
    257257        enddo ! iq=1,nq
    258258

  • trunk/LMDZ.PLUTO/libf/phypluto/convadj.F

    r3195 r3275  
    304304         end do
    305305         ! do l = 1, nlay
    306          !    print*,'zq(ig,:,vap)     = ',zq(i,l,igcm_h2o_vap)
     306         !    print*,'zq(ig,:,vap)     = ',zq(i,l,igcm_h2o_gas)
    307307         ! end do
    308308         ! do l = 1, nlay
    309          !    print*,'zq2(ig,:,vap)    = ',zq2(i,l,igcm_h2o_vap)
     309         !    print*,'zq2(ig,:,vap)    = ',zq2(i,l,igcm_h2o_gas)
    310310         ! end do
    311          !    print*,'zqm(vap)         = ',zqm(igcm_h2o_vap)
     311         !    print*,'zqm(vap)         = ',zqm(igcm_h2o_gas)
    312312            print*,'jadrs=',jadrs
    313313

  • trunk/LMDZ.PLUTO/libf/phypluto/dyn1d

    • Property svn:ignore set to
      *.mod

  • trunk/LMDZ.PLUTO/libf/phypluto/generic_cloud_common_h.F90

    r3184 r3275  
    1616    implicit none
    1717    real, save :: m  ! molecular mass of the specie (g/mol)
    18     real, save :: delta_vapH ! Enthalpy of vaporization (J/mol)
     18    real, save :: delta_gasH ! Enthalpy of vaporization (J/mol)
    1919    real,save :: Tref ! Ref temperature for Clausis-Clapeyron (K)
    2020    real,save :: Pref ! Reference pressure for Clausius Clapeyron (Pa)
    … …  
    4646        if (trim(specname) .eq. "Mg") then
    4747            print*,"Loading data for Mg"
    48             delta_vapH = 521.7E3 !J/mol
     48            delta_gasH = 521.7E3 !J/mol
    4949            Tref = 2303.0 !K
    5050            Pref =   1.0E5 !in Pa
    5151            m = 140.6931
    52             RLVTT_generic = delta_vapH/(m/1000.)
     52            RLVTT_generic = delta_gasH/(m/1000.)
    5353            metallicity_coeff=1.0*log(10.0)
    5454           
    5555        else if (trim(specname) .eq. "Na") then
    5656            print*,"Loading data for Na"
    57             delta_vapH = 265.9E3
     57            delta_gasH = 265.9E3
    5858            Tref = 1624.0
    5959            Pref =  1.0E5
    6060            m = 78.04454 !Na2S
    61             RLVTT_generic = delta_vapH/(m/1000.)
     61            RLVTT_generic = delta_gasH/(m/1000.)
    6262            metallicity_coeff=0.5*log(10.0)
    6363           
    6464        else if (trim(specname) .eq. "Fe") then
    6565            print*,"Loading data for Fe"
    66             delta_vapH = 401.9E3
     66            delta_gasH = 401.9E3
    6767            Tref = 2903.0
    6868            Pref  = 1.0E5
    6969            m = 55.8450
    70             RLVTT_generic = delta_vapH/(m/1000.)
     70            RLVTT_generic = delta_gasH/(m/1000.)
    7171            metallicity_coeff=0.0*log(10.0)
    7272           
    7373        else if (trim(specname) .eq. "Cr") then
    7474            print*,"Loading data for Cr"
    75             delta_vapH = 394.2E3
     75            delta_gasH = 394.2E3
    7676            Tref = 2749.0
    7777            Pref  = 1.0E5
    7878            m = 51.9961
    79             RLVTT_generic = delta_vapH/(m/1000.)
     79            RLVTT_generic = delta_gasH/(m/1000.)
    8080            metallicity_coeff=0.0*log(10.0)
    8181           
    8282        else if (trim(specname) .eq. "KCl") then
    8383            print*,"Loading data for KCl"
    84             delta_vapH = 217.9E3
     84            delta_gasH = 217.9E3
    8585            Tref = 1495.0
    8686            Pref = 1.0E5
    8787            metallicity_coeff=0.0*log(10.0)
    8888            m = 74.5498
    89             RLVTT_generic = delta_vapH/(m/1000.)
     89            RLVTT_generic = delta_gasH/(m/1000.)
    9090        else if (trim(specname) .eq. "Mn") then
    9191            print*,"Loading data for Mn"
    92             delta_vapH = 455.8E3
     92            delta_gasH = 455.8E3
    9393            Tref = 2064.0
    9494            Pref = 1.0E5
    9595            metallicity_coeff=1.0*log(10.0)
    9696            m = 87.003049
    97             RLVTT_generic = delta_vapH/(m/1000.)
     97            RLVTT_generic = delta_gasH/(m/1000.)
    9898        else if (trim(specname) .eq. "Zn") then
    9999            print*,"Loading data for Zn"
    100             delta_vapH = 303.9E3
     100            delta_gasH = 303.9E3
    101101            Tref = 1238.0
    102102            Pref = 1.0E5
    103103            metallicity_coeff=1.0*log(10.0)
    104104            m = 97.445
    105             RLVTT_generic = delta_vapH/(m/1000.)
     105            RLVTT_generic = delta_gasH/(m/1000.)
    106106        else
    107107            print*,"Unknow species (not in Mg, Fe, Na, KCl, Cr, Mn or Zn)"
    … …  
    138138            if (index(table_traceurs_line,specname) /= 0) then
    139139               
    140                 if (index(table_traceurs_line,'delta_vapH='   ) /= 0) then
    141                     read(table_traceurs_line(index(table_traceurs_line,'delta_vapH=')+len('delta_vapH='):),*) delta_vapH
    142                     print*, 'delta_vapH ', delta_vapH
     140                if (index(table_traceurs_line,'delta_gasH='   ) /= 0) then
     141                    read(table_traceurs_line(index(table_traceurs_line,'delta_gasH=')+len('delta_gasH='):),*) delta_gasH
     142                    print*, 'delta_gasH ', delta_gasH
    143143                end if
    144144                if (index(table_traceurs_line,'Tref='   ) /= 0) then
    … …  
    178178        close(117)
    179179
    180         RLVTT_generic=delta_vapH/(m/1000.)
     180        RLVTT_generic=delta_gasH/(m/1000.)
    181181
    182182        write(*,*) 'RLVTT_generic', RLVTT_generic
    … …  
    190190    !============================================================================
    191191    !   Clausius-Clapeyron relation :
    192     !   d(Ln(psat))/dT = delta_vapH/(RT^2)
    193     !   ->Psat = Pref * exp(-delta_vapH/R*(1/T-1/Tref))
     192    !   d(Ln(psat))/dT = delta_gasH/(RT^2)
     193    !   ->Psat = Pref * exp(-delta_gasH/R*(1/T-1/Tref))
    194194    !
    195195    !   Authors
    … …  
    203203        real, intent(out) :: psat,qsat !saturation vapor pressure (Pa) and mass mixing ratio at saturation (kg/kg) of the layer
    204204
    205         psat = pref*exp(-delta_vapH/(r*mugaz/1000.)*(1/T-1/Tref)-metallicity_coeff*metallicity) ! in Pa (because pref in Pa)
     205        psat = pref*exp(-delta_gasH/(r*mugaz/1000.)*(1/T-1/Tref)-metallicity_coeff*metallicity) ! in Pa (because pref in Pa)
    206206
    207207        if (psat .gt. p) then
    … …  
    217217    !===============================================================================
    218218    !   Compute dqsat = L/cp* d(qsat)/dT and d(Ln(psat)) = L/cp*d(Ln(psat))/dT
    219     !   we have d(ln(psat))/dT =  delta_vapH/R*(1/T^2) for clausius-clapeyron
     219    !   we have d(ln(psat))/dT =  delta_gasH/R*(1/T^2) for clausius-clapeyron
    220220    !   r*mugaz/1000. is the perfect gaz constant in the computation of "dummy"
    221221    !   Authors
    … …  
    235235        real :: dummy ! used to store d(ln(psat))/dT
    236236
    237         dummy = delta_vapH/((r*mugaz/1000.)*(T**2))
     237        dummy = delta_gasH/((r*mugaz/1000.)*(T**2))
    238238
    239239        if (psat .gt. p) then
    … …  
    270270
    271271
    272         Tsat = 1/(1/Tref - (r*mugaz/1000.)/delta_vapH*Log(p/Pref))
     272        Tsat = 1/(1/Tref - (r*mugaz/1000.)/delta_gasH*Log(p/Pref))
    273273
    274274    end subroutine Tsat_generic

  • trunk/LMDZ.PLUTO/libf/phypluto/generic_tracer_index_mod.F90

    r3184 r3275  
    66!     for both the vap & ice parts of the tracer
    77!
    8 !     call_ice_vap_generic is a boolean which tells if the tracer is condensable
    9 !     and allows, if it is condensable, to call the vap/ice pair 
     8!     call_ice_gas_generic is a boolean which tells if the tracer is condensable
     9!     and allows, if it is condensable, to call the vap/ice pair
    1010!
    1111!     Noé Clément & Lucas Teinturier (2022)
    … …  
    1717    contains
    1818
    19     subroutine generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
     19    subroutine generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
    2020
    2121        USE tracer_h, only: noms, is_condensable
    22         integer nq,iq, igcm_generic_vap, igcm_generic_ice, i_shortname, ii
    23         logical call_ice_vap_generic
     22        integer nq,iq, igcm_generic_gas, igcm_generic_ice, i_shortname, ii
     23        logical call_ice_gas_generic
    2424        character(len=10) :: shortname
    25         if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0) &
     25        if((is_condensable(iq)==1) .and. (index(noms(iq),"gas") .ne. 0) &
    2626                                                    .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
    27             ! Let's get the index of our tracers (we look for igcm _generic_vap and igcm_generic_ice)
     27            ! Let's get the index of our tracers (we look for igcm _generic_gas and igcm_generic_ice)
    2828
    29             igcm_generic_vap = iq
     29            igcm_generic_gas = iq
    3030            igcm_generic_ice = -1
    3131
    … …  
    4646            ! ! endif
    4747
    48             ! ! call_ice_vap_generic = .true.
     48            ! ! call_ice_gas_generic = .true.
    4949            i_shortname = index(noms(iq),'_')-1
    5050            shortname=noms(iq)(1:i_shortname)
    … …  
    5959            endif
    6060            ! if we didn't find it before, we look after the vap tracer in pq
    61             if (igcm_generic_ice .eq. -1 .and. (iq<nq)) then 
     61            if (igcm_generic_ice .eq. -1 .and. (iq<nq)) then
    6262                do ii=iq+1,nq
    63                     if (index(noms(ii),adjustl(trim(shortname))) .ne. 0) then 
     63                    if (index(noms(ii),adjustl(trim(shortname))) .ne. 0) then
    6464                        igcm_generic_ice=ii
    6565                        exit
    … …  
    6767                enddo
    6868            endif
    69             call_ice_vap_generic = .true.
     69            call_ice_gas_generic = .true.
    7070        else
    7171
    72             call_ice_vap_generic = .false.
     72            call_ice_gas_generic = .false.
    7373
    7474        end if ! if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0))

  • trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90

    r3258 r3275  
    883883      " albmin_ch4 = ",albmin_ch4
    884884
     885
     886     if (is_master)write(*,*)trim(rname)//&
     887     "fixed gaseous CH4 mixing ratio for rad transfer ?"
     888     ch4fix=.false. ! default value
     889     call getin_p("ch4fix",ch4fix)
     890     if (is_master)write(*,*)trim(rname)//&
     891       " ch4fix = ",ch4fix
     892     if (is_master)write(*,*)trim(rname)//&
     893       "fixed gaseous CH4 mixing ratio for rad transfer ?"
     894     vmrch4fix=0.5 ! default value
     895     call getin_p("vmrch4fix",vmrch4fix)
     896     if (is_master)write(*,*)trim(rname)//&
     897        " vmrch4fix = ",vmrch4fix
     898     vmrch4_proffix=.false. ! default value
     899     call getin_p("vmrch4_proffix",vmrch4_proffix)
     900     if (is_master)write(*,*)trim(rname)//&
     901        " vmrch4_proffix = ",vmrch4_proffix
     902
     903
    885904!*********** CO *********************************
    886905

  • trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90

    r3258 r3275  
    127127       IF (.NOT. allocated(is_rgcs)) allocate(is_rgcs(nq)) !LT
    128128       IF (.NOT. allocated(constants_mass)) allocate(constants_mass(nq))
    129        IF (.NOT. allocated(constants_delta_vapH)) allocate(constants_delta_vapH(nq))
     129       IF (.NOT. allocated(constants_delta_gasH)) allocate(constants_delta_gasH(nq))
    130130       IF (.NOT. allocated(constants_Tref)) allocate(constants_Tref(nq))
    131131       IF (.NOT. allocated(constants_Pref)) allocate(constants_Pref(nq))
    … …  
    157157       is_rgcs(:) = 0
    158158       constants_mass(:)=0
    159        constants_delta_vapH(:)=0
     159       constants_delta_gasH(:)=0
    160160       constants_Tref(:)=0
    161161       constants_Pref(:)=0
    … …  
    229229          write(*,*) 'Tracer ',count,' = n2'
    230230        endif
    231         if (noms(iq).eq."ch4_vap") then
     231        if (noms(iq).eq."ch4_gas") then
    232232          igcm_ch4_gas=iq
    233233          mmol(igcm_ch4_gas)=16.
    … …  
    243243          write(*,*) 'Tracer ',count,' = ch4 ice'
    244244        endif
    245         if (noms(iq).eq."co_vap") then
     245        if (noms(iq).eq."co_gas") then
    246246          igcm_co_gas=iq
    247247          mmol(igcm_co_gas)=28.
    … …  
    376376                call specie_parameters_table(noms(iq)(1:len(trim(noms(iq)))-4))
    377377                constants_mass(iq)=m
    378                 constants_delta_vapH(iq)=delta_vapH
     378                constants_delta_gasH(iq)=delta_gasH
    379379                constants_Tref(iq)=Tref
    380380                constants_Pref(iq)=Pref

  • trunk/LMDZ.PLUTO/libf/phypluto/mass_redistribution_mod.F90

    r3184 r3275  
    77                       pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr,  &
    88                       pdtmr,pdtsrfmr,pdpsrfmr,pdumr,pdvmr,pdqmr,pdqsmr)
    9                                                    
     9
    1010       use radcommon_h, only: glat
    1111       USE planete_mod, only: bp
    1212       use comcstfi_mod, only: g
    13        
     13
    1414       IMPLICIT NONE
    1515!=======================================================================
    … …  
    2222!
    2323!   input:
    24 !   ----- 
     24!   -----
    2525!    ngrid                 nombre de points de verticales
    2626!                          (toutes les boucles de la physique sont au
    2727!                          moins vectorisees sur ngrid)
    2828!    nlayer                nombre de couches
    29 !    pplay(ngrid,nlayer)   Pressure levels 
    30 !    pplev(ngrid,nlayer+1) Pressure levels 
     29!    pplay(ngrid,nlayer)   Pressure levels
     30!    pplev(ngrid,nlayer+1) Pressure levels
    3131!    nq                    Number of tracers
    3232!
    … …  
    3535!    pu,pv (ngrid,nlayer)  wind velocity (m/s)
    3636!
    37 !                   
     37!
    3838!    pdX                   physical tendency of X before mass redistribution
    3939!
    … …  
    4545!    pdXmr(ngrid)           physical tendency of X after mass redistribution
    4646!
    47 !   
     47!
    4848!
    4949!=======================================================================
    … …  
    5555!    Arguments :
    5656!    ---------
    57       INTEGER,INTENT(IN) :: ngrid, nlayer, nq   
     57      INTEGER,INTENT(IN) :: ngrid, nlayer, nq
    5858      REAL,INTENT(IN) :: ptimestep
    5959      REAL,INTENT(IN) :: pcapcal(ngrid)
    … …  
    8080      REAL zzu(nlayer),zzv(nlayer)
    8181      REAL zzq(nlayer,nq),zzt(nlayer)
    82 !    Dummy variables     
     82!    Dummy variables
    8383      INTEGER n,l,ig,iq
    8484      REAL zdtsig(ngrid,nlayer)
    … …  
    8888      REAL zq1(nlayer)
    8989      REAL ztsrf(ngrid)
    90       REAL ztm(nlayer+1) 
     90      REAL ztm(nlayer+1)
    9191      REAL zum(nlayer+1) , zvm(nlayer+1)
    9292      REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
    … …  
    103103!
    104104      IF (firstcall) THEN
    105          firstcall=.false.       
     105         firstcall=.false.
    106106      ENDIF
    107107!
    108108!======================================================================
    109 !    Calcul of h2o condensation 
     109!    Calcul of h2o condensation
    110110!    ============================================================
    111 ! 
     111!
    112112!    Used variable :
    113113!       pdmassmr      : air Mass added to the atmosphere in each layer per unit time (kg/m2/s)
    … …  
    135135!    Changing pressure at the surface:
    136136!    """"""""""""""""""""""""""""""""""""
    137          
     137
     138      zmassboil(1:ngrid) = 0 ! AF: TODO: update this for pluto
     139      pdtsrfmr(1:ngrid) = 0
     140
    138141      pdpsrfmr(1:ngrid) = (zdmass_sum(1:ngrid,1)+zmassboil(1:ngrid))*g
    139142
    … …  
    151154
    152155! ***************************************************************
    153 !  Correction to account for redistribution between sigma or hybrid 
     156!  Correction to account for redistribution between sigma or hybrid
    154157!  layers when changing surface pressure
    155158!  zzx quantities have dimension (nlayer) to speed up calculation
    … …  
    174177           ! Ehouarn: but for now leave things as before
    175178            zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - (bp(l)-bp(l+1))*(zdmass_sum(ig,1)+zmflux(1))
    176 ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld 
     179! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
    177180            if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
    178181         END DO
    179182
    180183! Mass of each layer
    181 ! ------------------ 
     184! ------------------
    182185         zzmass(1:nlayer)=zmass(ig,1:nlayer)*(1.+pdpsrfmr(ig)*ptimestep/pplev(ig,1))
    183186
    … …  
    186189!  """"""""""""""""""""""""""""""""
    187190
    188 !        averaging operator for TRANSPORT 
     191!        averaging operator for TRANSPORT
    189192!        """"""""""""""""""""""""""""""""
    190193!        Value transfert at the surface interface when condensation
    191194!        sublimation:
    192195         ztm(1) = ztsrf(ig)
    193          zum(1) = 0. 
    194          zvm(1) = 0. 
     196         zum(1) = 0.
     197         zvm(1) = 0.
    195198         zqm(1,1:nq)=0. ! most tracer do not condense !
    196          !! if (water) zqm(1,igcm_h2o_vap)=1. ! flux is 100% h2o at surface
    197          
     199         !! if (water) zqm(1,igcm_h2o_gas)=1. ! flux is 100% h2o at surface
     200
    198201!        Van Leer scheme:
    199202         w(1:nlayer+1)=-zmflux(1:nlayer+1)*ptimestep
    200          call vl1d(nlayer,zzt,2.,zzmass,w,ztm) 
    201          call vl1d(nlayer,zzu,2.,zzmass,w,zum) 
    202          call vl1d(nlayer,zzv,2.,zzmass,w,zvm) 
     203         call vl1d(nlayer,zzt,2.,zzmass,w,ztm)
     204         call vl1d(nlayer,zzu,2.,zzmass,w,zum)
     205         call vl1d(nlayer,zzv,2.,zzmass,w,zvm)
    203206         do iq=1,nq
    204207           zq1(1:nlayer)=zzq(1:nlayer,iq)
    … …  
    214217
    215218!        Surface condensation affects low winds
    216          if (zmflux(1).lt.0) then 
     219         if (zmflux(1).lt.0) then
    217220            zum(1)= zzu(1) *  (w(1)/zzmass(1))
    218221            zvm(1)= zzv(1) *  (w(1)/zzmass(1))
    … …  
    222225            end if
    223226         end if
    224                    
    225          ztm(nlayer+1)= zzt(nlayer) ! should not be used, but... 
     227
     228         ztm(nlayer+1)= zzt(nlayer) ! should not be used, but...
    226229         zum(nlayer+1)= zzu(nlayer)  ! should not be used, but...
    227230         zvm(nlayer+1)= zzv(nlayer)  ! should not be used, but...
    228231         zqm(nlayer+1,1:nq)= zzq(nlayer,1:nq)
    229  
    230 !        Tendencies on T, U, V, Q 
     232
     233!        Tendencies on T, U, V, Q
    231234!        """"""""""""""""""""""""
    232235         DO l=1,nlayer
    233  
     236
    234237!           Tendencies on T
    235238            pdtmr(ig,l) = (1/zzmass(l)) *   &
    … …  
    253256         enddo
    254257
    255       END DO  ! loop on ig 
     258      END DO  ! loop on ig
    256259
    257260CONTAINS
    … …  
    260263      SUBROUTINE vl1d(llm,q,pente_max,zzmass,w,qm)
    261264!
    262 !   
     265!
    263266!     Operateur de moyenne inter-couche pour calcul de transport type
    264267!     Van-Leer " pseudo amont " dans la verticale
    … …  
    270273!
    271274!   --------------------------------------------------------------------
    272      
     275
    273276      IMPLICIT NONE
    274277
    … …  
    280283      REAL w(llm+1)
    281284!
    282 !      Local 
     285!      Local
    283286!   ---------
    284287!
    … …  
    287290      real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
    288291      real sigw, Mtot, MQtot
    289       integer m 
    290 !     integer ismax,ismin 
    291 
    292 
    293 !    On oriente tout dans le sens de la pression 
     292      integer m
     293!     integer ismax,ismin
     294
     295
     296!    On oriente tout dans le sens de la pression
    294297!     W > 0 WHEN DOWN !!!!!!!!!!!!!
    295298
    … …  
    347350               stop
    348351             end if
    349           else      ! if(w(l+1).lt.0) 
    350              m = l-1 
     352          else      ! if(w(l+1).lt.0)
     353             m = l-1
    351354             Mtot = zzmass(m+1)
    352355             MQtot = zzmass(m+1)*q(m+1)
    … …  
    365368             else
    366369                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
    367              end if   
     370             end if
    368371          end if
    369372       enddo
    … …  
    373376!         if(w(1).gt.0.) then
    374377!            qm(1)=q(1)
    375 !         else 
     378!         else
    376379!           qm(1)=0.
    377380!         end if

  • trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90

    r3184 r3275  
    257257endif ! of if (nq.ge.1)
    258258
    259 !!call WriteField_phy("in_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     259!!call WriteField_phy("in_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
    260260
    261261if (startphy_file) then

  • trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

    r3258 r3275  
    300300      real omega(ngrid,nlayer)            ! omega velocity (Pa/s, >0 when downward)
    301301
    302       integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice
    303       logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
     302      integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_gas, igcm_generic_ice
     303      logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer
    304304
    305305      real zls                       ! Solar longitude (radians).
    … …  
    638638         qsurf_hist(:,:)=qsurf(:,:)
    639639
    640 !!         call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     640!!         call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
    641641
    642642!        Initialize variable for dynamical heating and zonal wind tendency diagnostic
    … …  
    681681         endif
    682682
    683 !!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     683!!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
    684684#endif
    685685         if (corrk) then
    … …  
    702702         endif
    703703
    704 !!         call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     704!!         call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
    705705         ! XIOS outputs
    706706#ifdef CPP_XIOS
    … …  
    711711#endif
    712712
    713 !!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     713!!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
    714714
    715715         write(*,*) "physiq: end of firstcall"
    716716      endif ! end of 'firstcall'
    717717
    718 !!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    719 !!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     718!!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
     719!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
    720720
    721721      if (check_physics_inputs) then
    … …  
    882882            ! Eclipse incoming sunlight !AF24: removed
    883883
    884 !!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
    885 !!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     884!!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
     885!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
    886886
    887887
    … …  
    900900                  do iq=1,nq
    901901
    902                      call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    903 
    904                      if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
     902                     call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
     903
     904                     if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
    905905
    906906                        epsi_generic=constants_epsi_generic(iq)
    907907
    908                         muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
    909                         muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
     908                        muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_gas))
     909                        muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_gas))
    910910
    911911                     endif
    … …  
    931931                call callcorrk(ngrid,nlayer,pq,nq,qsurf,                          &
    932932                              albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
    933                               tsurf,fract,dist_star,aerosol,muvar,                &
     933                              zzlay,tsurf,fract,dist_star,aerosol,muvar,                &
    934934                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
    935935                              fluxsurfabs_sw,fluxtop_lw,                          &
    … …  
    10391039      endif ! of if (callrad)
    10401040
    1041 !!      call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
    1042 !!      call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     1041!!      call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
     1042!!      call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
    10431043
    10441044
    … …  
    11371137         ! endif
    11381138
    1139 !!         call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_vap))
     1139!!         call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_gas))
    11401140
    11411141         if (tracer) then
    … …  
    11441144         end if ! of if (tracer)
    11451145
    1146 !!         call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
    1147 !!         call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     1146!!         call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
     1147!!         call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
    11481148
    11491149         ! test energy conservation
    … …  
    12641264         dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2)
    12651265
    1266 !!         call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
    1267 !!         call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     1266!!         call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
     1267!!         call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
    12681268
    12691269         ! test energy conservation
    … …  
    14581458               cloudfrac(:,:)=0.0
    14591459               do iq=1,nq
    1460                   call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    1461                   if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
     1460                  call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
     1461                  if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
    14621462                     do l = 1, nlayer
    14631463                        do ig=1,ngrid
    … …  
    14911491            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
    14921492
    1493 !!            call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1493!!            call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
    14941494
    14951495            ! ! Test water conservation !AF24: removed
    … …  
    14971497         end if ! end of 'sedimentation'
    14981498
    1499 !!         call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
    1500 !!         call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     1499!!         call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
     1500!!         call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
    15011501
    15021502  ! ---------------
    … …  
    15091509            zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0
    15101510            !    (   zdqevap(1:ngrid,1:nlayer)                          &
    1511             !    !   + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)             &
    1512             !    !   + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)             &
     1511            !    !   + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas)             &
     1512            !    !   + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas)             &
    15131513            !      + dqvaplscale(1:ngrid,1:nlayer) )
    15141514
    … …  
    15361536         endif
    15371537
    1538 !         call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
    1539 
    1540 !!         call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
    1541 !!         call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     1538!         call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
     1539
     1540!!         call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
     1541!!         call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
    15421542
    15431543
    … …  
    17941794         do iq=1,nq
    17951795
    1796             call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    1797 
    1798             if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
     1796            call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
     1797
     1798            if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
    17991799
    18001800               do l = 1, nlayer
    18011801                  do ig=1,ngrid
    18021802                     call Psat_generic(zt(ig,l),pplay(ig,l),metallicity,psat_tmp_generic,qsat_generic(ig,l,iq))
    1803                      RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_vap) / qsat_generic(ig,l,iq)
     1803                     RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_gas) / qsat_generic(ig,l,iq)
    18041804                  enddo
    18051805               enddo
    … …  
    23642364
    23652365      do iq=1,nq
    2366          call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    2367 
    2368          if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
     2366         call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
     2367
     2368         if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
    23692369
    23702370            reffrad_generic_zeros_for_wrf(:,:) = 1.
    … …  
    23732373            comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
    23742374            comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
    2375             comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_vap)
     2375            comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_gas)
    23762376            comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
    23772377            ! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
    … …  
    24172417!        do ig=1,ngrid
    24182418!          if ((tracer).and.(water)) then
    2419 !           pdq(ig,:,igcm_h2o_vap) = pdq(ig,:,igcm_h2o_vap) + lsf_dq(:)
     2419!           pdq(ig,:,igcm_h2o_gas) = pdq(ig,:,igcm_h2o_gas) + lsf_dq(:)
    24202420!          endif
    24212421!          pdt(ig,:) = pdt(ig,:) + lsf_dt(:)

  • trunk/LMDZ.PLUTO/libf/phypluto/sfluxi.F

    r3184 r3275  
    11      module sfluxi_mod
    2      
     2
    33      implicit none
    4      
     4
    55      contains
    6      
     6
    77      SUBROUTINE SFLUXI(PLEV,TLEV,DTAUI,TAUCUMI,UBARI,RSFI,WNOI,DWNI,
    88     *                  COSBI,WBARI,NFLUXTOPI,NFLUXTOPI_nu,
    9      *                  FMNETI,fluxupi,fluxdni,fluxupi_nu,
     9     *                  FMNETI,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,
    1010     *                  FZEROI,TAUGSURF)
    11      
     11
    1212      use radinc_h, only: NTfac, NTstart, L_LEVELS, L_NSPECTI, L_NGAUSS
    1313      use radinc_h, only: L_NLAYRAD, L_NLEVRAD
    … …  
    1515      use comcstfi_mod, only: pi
    1616      use gfluxi_mod, only: gfluxi
    17      
     17
    1818      implicit none
    19      
     19
    2020      integer NLEVRAD, L, NW, NG, NTS, NTT
    21      
     21
    2222      real*8 TLEV(L_LEVELS), PLEV(L_LEVELS)
    2323      real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
    2424      real*8 FMNETI(L_NLAYRAD)
     25      real*8 FMNETI_NU(L_NLAYRAD,L_NSPECTI)
    2526      real*8 WNOI(L_NSPECTI), DWNI(L_NSPECTI)
    2627      real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    … …  
    3233      real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)
    3334      real*8 FTOPUP
    34      
     35
    3536      real*8 UBARI, RSFI, TSURF, BSURF, TTOP, BTOP, TAUTOP
    3637      real*8 PLANCK, PLTOP
    … …  
    3839      real*8 FZEROI(L_NSPECTI)
    3940      real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero
    40      
     41
    4142      real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
    4243      real*8 PLANCKSUM,PLANCKREF
    43      
     44
    4445! AB : variables for interpolation
    4546      REAL*8 C1
    4647      REAL*8 C2
    4748      REAL*8 P1
    48      
     49
    4950!======================================================================C
    50      
     51
    5152      NLEVRAD = L_NLEVRAD
    52      
     53
    5354! ZERO THE NET FLUXES
    5455      NFLUXTOPI = 0.0D0
    55      
     56
    5657      DO NW=1,L_NSPECTI
    5758        NFLUXTOPI_nu(NW) = 0.0D0
    5859        DO L=1,L_NLAYRAD
    5960           FLUXUPI_nu(L,NW) = 0.0D0
     61           FMNETI_nu(L,NW)  = 0.0
    6062           fup_tmp(nw)=0.0D0
    6163           fdn_tmp(nw)=0.0D0
    6264        END DO
    6365      END DO
    64      
     66
    6567      DO L=1,L_NLAYRAD
    6668        FMNETI(L)  = 0.0D0
    … …  
    6870        FLUXDNI(L) = 0.0D0
    6971      END DO
    70      
     72
    7173! WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED
    7274! TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
    73      
     75
    7476      TTOP  = TLEV(2)  ! JL12 why not (1) ???
    7577      TSURF = TLEV(L_LEVELS)
    … …  
    8688      PLANCKREF = TSURF * TSURF
    8789      PLANCKREF = sigma * PLANCKREF * PLANCKREF
    88      
     90
    8991      DO NW=1,L_NSPECTI
    9092! AB : PLANCKIR(NW,NTS) is replaced by P1, the linear interpolation result for a temperature TSURF
    … …  
    9395         PLANCKSUM = PLANCKSUM + P1 * DWNI(NW)
    9496      ENDDO
    95      
     97
    9698      PLANCKSUM = PLANCKREF / (PLANCKSUM * Pi)
    9799!JL12
    98      
     100
    99101      DO 501 NW=1,L_NSPECTI
    100102! SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
    … …  
    106108         BSURF = (1. - RSFI) * P1 * PLANCKSUM
    107109         PLTOP = (1.0D0 - C2) * PLANCKIR(NW,NTT) + C2*PLANCKIR(NW,NTT+1)
    108          
     110
    109111! If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the
    110112! special Gauss point at the end.
    111113         FZERO = FZEROI(NW)
    112          
     114
    113115         IF(FZERO.ge.0.99) goto 40
    114          
     116
    115117         DO NG=1,L_NGAUSS-1
    116            
     118
    117119            if(TAUGSURF(NW,NG).lt. TLIMIT) then
    118120               fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG)
    119121               goto 30
    120122            end if
    121            
     123
    122124! SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR
    123125! CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
    124126! OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
    125            
     127
    126128!            TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
    127129            TAUTOP = TAUCUMI(2,NW,NG)
    128130            BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
    129            
     131
    130132! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
    131133! CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
    132 ! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 
    133            
     134! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
     135
    134136            CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
    135137     *                TAUCUMI(1,NW,NG),
    136138     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
    137139     *                BSURF,FTOPUP,FMUPI,FMDI)
    138          
     140
    139141! NOW CALCULATE THE CUMULATIVE IR NET FLUX
    140142            NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)
    141143     *                * (1.0D0-FZEROI(NW))
    142            
     144
    143145! and same thing by spectral band... (RDW)
    144146            NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW) + FTOPUP * DWNI(NW)
    145147     *                       * GWEIGHT(NG) * (1.0D0-FZEROI(NW))
    146            
     148
    147149            DO L=1,L_NLEVRAD-1
    148150!           CORRECT FOR THE WAVENUMBER INTERVALS
    149151               FMNETI(L)  = FMNETI(L) + (FMUPI(L)-FMDI(L)) * DWNI(NW)
    150152     *                    * GWEIGHT(NG)*(1.0D0-FZEROI(NW))
     153               FMNETI_NU(L,NW)  = FMNETI_NU(L,NW) + (FMUPI(L)-FMDI(L))
     154     *                    * DWNI(NW) * GWEIGHT(NG)*(1.0D0-FZEROI(NW))
    151155               FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG)
    152156     *                    * (1.0D0-FZEROI(NW))
    … …  
    157161     *                          * GWEIGHT(NG) * (1.0D0 - FZEROI(NW))
    158162            END DO
    159            
     163
    160164   30       CONTINUE
    161          
     165
    162166         END DO       !End NGAUSS LOOP
    163          
     167
    164168   40    CONTINUE
    165          
     169
    166170! SPECIAL 17th Gauss point
    167171         NG     = L_NGAUSS
    168          
     172
    169173!         TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
    170174         TAUTOP = TAUCUMI(2,NW,NG)
    171175         BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
    172          
     176
    173177! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
    174178! CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
    175 ! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 
    176          
     179! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
     180
    177181         CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
    178182     *                TAUCUMI(1,NW,NG),
    179183     *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
    180184     *                BSURF,FTOPUP,FMUPI,FMDI)
    181          
     185
    182186! NOW CALCULATE THE CUMULATIVE IR NET FLUX
    183187         NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*FZERO
    184          
     188
    185189!         and same thing by spectral band... (RW)
    186190         NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
    187191     *      +FTOPUP*DWNI(NW)*FZERO
    188          
     192
    189193         DO L=1,L_NLEVRAD-1
    190194! CORRECT FOR THE WAVENUMBER INTERVALS
    191195            FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*FZERO
     196            FMNETI_nu(L,NW)  = FMNETI_nu(L,NW)+(FMUPI(L)-FMDI(L))
     197     *       *DWNI(NW)*FZERO
    192198            FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*FZERO
    193199            FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*FZERO
    … …  
    196202     *                       + FMUPI(L) * DWNI(NW) * FZERO
    197203         END DO
    198          
     204
    199205  501 CONTINUE      !End Spectral Interval LOOP
    200206! *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****
    201      
     207
    202208      END SUBROUTINE SFLUXI
    203209

  • trunk/LMDZ.PLUTO/libf/phypluto/sfluxv.F

    r3184 r3275  
    22
    33      implicit none
    4      
     4
    55      contains
    6      
     6
    77      SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV,
    88     *                  UBAR0,STEL,NFLUXTOPV,FLUXTOPVDN,
    9      *                  NFLUXOUTV_nu,NFLUXGNDV_nu,
    10      *                  FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf)
     9     *                  NFLUXOUTV_nu,NFLUXGNDV_nu,NFLUXTOPV_nu,
     10     *                  FMNETV,FMNETV_NU,FLUXUPV,FLUXDNV,
     11     *                  FZEROV,taugsurf)
    1112
    1213      use radinc_h, only: L_TAUMAX, L_LEVELS, L_NSPECTV, L_NGAUSS
    … …  
    1819
    1920      real*8 FMNETV(L_NLAYRAD)
     21      real*8 FMNETV_NU(L_NLAYRAD,L_NSPECTV)
    2022      real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
    2123      real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
    … …  
    2729      real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD)
    2830      real*8 NFLUXTOPV, FLUXUP, FLUXDN,FLUXTOPVDN
     31      real*8 NFLUXTOPV_nu(L_NSPECTV)
    2932      real*8 NFLUXOUTV_nu(L_NSPECTV)
    3033      real*8 NFLUXGNDV_nu(L_NSPECTV)
    … …  
    5053         NFLUXOUTV_nu(NW)=0.0
    5154         NFLUXGNDV_nu(NW)=0.0
     55         NFLUXTOPV_nu(nw)=0.0
     56         DO L=1,L_NLAYRAD
     57           FMNETV_NU(L,NW) = 0.0
     58        END DO
    5259      END DO
    5360
    … …  
    6471
    6572      DO 500 NW=1,L_NSPECTV
    66      
     73
    6774        F0PI = STEL(NW)
    6875
    … …  
    8491          BSURF = 0.
    8592C         LOOP OVER THE NTERMS BEGINNING HERE
    86  
     93
    8794
    8895!      FACTOR    = 1.0D0 - WDEL(1)*CDEL(1)**2
    … …  
    96103C         CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
    97104C         WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
    98 C 
    99 C         FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO 
     105C
     106C         FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO
    100107C         RETURN FLUXES FOR A GIVEN NT
    101108
    102109
    103110          CALL GFLUXV(DTAUV(1,NW,NG),TAUV(1,NW,NG),TAUCUMV(1,NW,NG),
    104      *                WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW),   
     111     *                WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW),
    105112     *                BTOP,BSURF,FMUPV,FMDV,DIFFV,FLUXUP,FLUXDN)
    106113
    107 C         NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX 
     114C         NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX
    108115
    109116          NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*GWEIGHT(NG)*
    … …  
    114121            FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )*
    115122     *                           GWEIGHT(NG)*(1.0-FZEROV(NW))
     123            FMNETV_NU(L,NW)=FMNETV_NU(L,NW)+( FMUPV(L)-FMDV(L) )*
     124     *                           GWEIGHT(NG)*(1.0-FZEROV(NW))
    116125            FLUXUPV(L) = FLUXUPV(L) + FMUPV(L)*GWEIGHT(NG)*
    117126     *                   (1.0-FZEROV(NW))
    … …  
    120129          END DO
    121130
     131c         and same thing by spectral band... (RDW)
     132          NFLUXTOPV_nu(NW) = NFLUXTOPV_nu(NW)
     133     *      +(FLUXUP-FLUXDN)*GWEIGHT(NG)*
     134     *                          (1.0-FZEROV(NW))
     135
    122136c     band-resolved flux leaving TOA (RDW)
    123137          NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
    … …  
    133147          DIFFVT = DIFFVT + DIFFV*GWEIGHT(NG)*(1.0-FZEROV(NW))
    134148
    135    30     CONTINUE 
    136 
    137         END DO   ! the Gauss loop 
    138 
    139    40   continue 
     149   30     CONTINUE
     150
     151        END DO   ! the Gauss loop
     152
     153   40   continue
    140154C       Special 17th Gauss point
    141155
    … …  
    143157
    144158C       SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE VISIBLE
    145  
     159
    146160        BTOP = 0.0
    147161
    148162C       LOOP OVER THE NTERMS BEGINNING HERE
    149  
     163
    150164        ETERM = MIN(TAUV(L_NLEVRAD,NW,NG),TAUMAX)
    151165        BSURF = RSFV(NW)*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0)
    … …  
    155169C       CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
    156170C       WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
    157 C 
    158 C       FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO 
     171C
     172C       FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO
    159173C       RETURN FLUXES FOR A GIVEN NT
    160174
    … …  
    164178
    165179
    166 C       NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX 
     180C       NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX
    167181
    168182        NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*FZERO
    … …  
    170184        DO L=1,L_NLAYRAD
    171185          FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )*FZERO
     186          FMNETV_NU(L,NW)=FMNETV_NU(L,NW)+( FMUPV(L)-FMDV(L) )*FZERO
    172187          FLUXUPV(L) = FLUXUPV(L) + FMUPV(L)*FZERO
    173188          FLUXDNV(L) = FLUXDNV(L) + FMDV(L)*FZERO
    174189        END DO
    175190
     191c         and same thing by spectral band... (RDW)
     192          NFLUXTOPV_nu(NW) = NFLUXTOPV_nu(NW)
     193     *      +(FLUXUP-FLUXDN)*FZERO
     194
    176195c     band-resolved flux leaving TOA (RDW)
    177196          NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
    … …  
    196215
    197216      end module sfluxv_mod
    198      
     217

  • trunk/LMDZ.PLUTO/libf/phypluto/surfdat_h.F90

    r3195 r3275  
    1717!$OMP THREADPRIVATE(zmea,zstd,zsig,zgam,zthe)
    1818       real ttop
    19        real,allocatable,dimension(:) :: kp ! TB ref pressure 
     19       real,allocatable,dimension(:) :: kp ! TB ref pressure
    2020       real p00
    2121!$OMP THREADPRIVATE(ttop,kp,p00)
    2222! surface properties ! TB16
    2323       real alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a
    24 !$OMP THREADPRIVATE(alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a)       
     24!$OMP THREADPRIVATE(alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a)
    2525       real emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq
    26 !$OMP THREADPRIVATE(emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq)       
     26!$OMP THREADPRIVATE(emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq)
    2727       real ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe
    28 !$OMP THREADPRIVATE(ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe)       
     28!$OMP THREADPRIVATE(ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe)
    2929       real alb_tho_spe,emis_tho_spe
    30 !$OMP THREADPRIVATE(alb_tho_spe,emis_tho_spe)       
     30!$OMP THREADPRIVATE(alb_tho_spe,emis_tho_spe)
    3131
    3232       real,allocatable,dimension(:) :: dryness  !"Dryness coefficient" for grnd water ice sublimation

  • trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90

    r3237 r3275  
    5151       ! Lists of constants for condensable tracers
    5252       real, save, allocatable :: constants_mass(:)                    ! molecular mass of the specie (g/mol)
    53        real, save, allocatable :: constants_delta_vapH(:)           ! Enthalpy of vaporization (J/mol)
     53       real, save, allocatable :: constants_delta_gasH(:)           ! Enthalpy of vaporization (J/mol)
    5454       real, save, allocatable :: constants_Tref(:)                 ! Ref temperature for Clausis-Clapeyron (K)
    5555       real, save, allocatable :: constants_Pref(:)                 ! Reference pressure for Clausius Clapeyron (Pa)
    … …  
    5858       real, save, allocatable :: constants_metallicity_coeff(:)    ! Coefficient to take into account the metallicity
    5959       real, save, allocatable :: constants_RCPV_generic(:)                   ! specific heat capacity of the tracer vapor at Tref
    60 !$OMP THREADPRIVATE(constants_mass,constants_delta_vapH,constants_Tref)
     60!$OMP THREADPRIVATE(constants_mass,constants_delta_gasH,constants_Tref)
    6161!$OMP THREADPRIVATE(constants_Pref,constants_epsi_generic)
    6262!$OMP THREADPRIVATE(constants_RLVTT_generic,constants_metallicity_coeff,constants_RCPV_generic)

  • trunk/LMDZ.PLUTO/libf/phypluto/vdifc_mod.F

    r3195 r3275  
    435435         do iq=1,nq
    436436
    437             ! if (iq.ne.igcm_h2o_vap) then
     437            ! if (iq.ne.igcm_h2o_gas) then
    438438
    439439               DO ig=1,ngrid
Note: See TracChangeset for help on using the changeset viewer.