Ignore:
Timestamp:
May 31, 2024, 10:03:40 PM (6 months ago)
Author:
afalco
Message:

Pluto PCM:
Added zrecast & old sedim ;
Choose haze file ;
AF

Location:
trunk/LMDZ.PLUTO/libf/phypluto
Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

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

    r3196 r3353  
    178178      INTEGER :: ngrid,nlayer
    179179!     Aerosol effective radius used for radiative transfer (meter)
    180       REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind)
     180      REAL :: reffrad(ngrid,nlayer,naerkind)
    181181!     Aerosol effective variance used for radiative transfer (n.u.)
    182182      REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind)
    … …  
    324324        minrad=min(MINVAL(radiustab(1,1,1:nsize(1,1))),MINVAL(radiustab(1,2,1:nsize(1,2))))
    325325        maxrad=min(MAXVAL(radiustab(1,1,1:nsize(1,1))),MAXVAL(radiustab(1,2,1:nsize(1,2))))
    326         IF ((MINVAL(reffrad).LE.minrad).OR.(MAXVAL(reffrad).GE.maxrad)) then
     326        IF ((MINVAL(reffrad).LT.minrad).OR.(MAXVAL(reffrad).GT.maxrad)) then
    327327           WRITE(*,*) 'Warning: particle size in grid box #'
    328            WRITE(*,*) ig,' is too large to be used by the '
    329            WRITE(*,*) 'radiative transfer; please extend the '
    330            WRITE(*,*) 'interpolation grid to larger grain sizes.'
     328           WRITE(*,*) ig,' is larger than optical properties. '
     329           WRITE(*,*) 'reffrad=',MINVAL(reffrad),'-',MAXVAL(reffrad)
    331330           WRITE(*,*) 'radiustab=',minrad,'-',maxrad
    332            WRITE(*,*) 'reffrad=',MINVAL(reffrad),'-',MAXVAL(reffrad)
    333            stop
     331
     332           ! ensure reffrad is within bounds of radiustab
     333           WHERE(reffrad.LT.minrad)
     334              reffrad=minrad
     335           ENDWHERE
     336           WHERE(reffrad.GT.maxrad)
     337              reffrad=maxrad
     338           ENDWHERE
     339           WRITE(*,*) 'Truncated reffrad within radiustab bounds:'
     340           WRITE(*,*) 'New reffrad=',MINVAL(reffrad),'-',MAXVAL(reffrad)
    334341        ENDIF
    335342

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

    r3329 r3353  
    6767      parameter(Nfine=701)
    6868      character(len=100) :: file_path
     69      character(len=100) :: file_name
    6970      real,save :: levdat(Nfine),densdat(Nfine)
    7071
    … …  
    7475      IF (firstcall) then
    7576        firstcall=.false.
    76         file_path=trim(datadir)//'/haze_prop/'//hazemmr_file
     77
     78        if (hazemmr_file/='None')then
     79            file_name = hazemmr_file
     80            print*, 'Read Haze MMR file: ',hazemmr_file
     81        else if(hazedens_file/='None')then
     82            file_name = hazedens_file
     83            print*, 'Read Haze density: ',hazedens_file
     84        else
     85            STOP "No filename given for haze profile. Either set hazemmr_file or hazedens_file"
     86        endif
     87
     88        file_path=trim(datadir)//'/haze_prop/'//file_name
    7789        open(224,file=file_path,form='formatted')
    7890        do ifine=1,Nfine
    … …  
    8092        enddo
    8193        close(224)
    82         print*, 'TB22 read Haze MMR profile'
     94        print*, 'Read Haze profile: ',file_path
    8395      ENDIF
    8496
    … …  
    91103      !                                --> kg m-3 --> kg/kg
    92104      do iaer=1,naerkind
    93             if(iaer.eq.iaero_haze.and.1.eq.2) then !TB22 activate/deactivate mmr or part density
     105            if(iaer.eq.iaero_haze.and.hazedens_file/='None') then !AF24 activate/deactivate mmr or part density
    94106              !print*, 'Haze profile is fixed'
    95107              do ig=1,ngrid
    … …  
    97109                    !from number density in cm-3
    98110                    profmmr(ig,l)=profmmr(ig,l)*1.e6*4./3.*pi*reffrad(ig,l,iaer)**3*rho_q(i_haze)/(pplay(ig,l)/(r*pt(ig,l)))
     111                ! print*, profmmr(ig,l)
    99112                 enddo
    100113              enddo

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

    r3184 r3353  
    1111!
    1212!     Authors
    13 !     ------- 
     13!     -------
    1414!     Robin Wordsworth (2009)
    1515!     A. Spiga: make the routine OK with latest changes in rcm1d
    … …  
    2525
    2626      character(len=20),parameter :: myname="calc_cpp_mugaz"
    27       real cpp_c   
     27      real cpp_c
    2828      real mugaz_c
    2929
    … …  
    4141         else
    4242            ! all values at 300 K from Engineering Toolbox
    43             if(igas.eq.igas_N2)then
     43            if(igas.eq.igas_CO2)then
    4444               mugaz_c = mugaz_c + 44.01*gfrac(igas)
    4545            elseif(igas.eq.igas_N2)then
    … …  
    5959            elseif(igas.eq.igas_NH3)then
    6060               mugaz_c = mugaz_c + 17.03*gfrac(igas)
    61             elseif(igas.eq.igas_C2H6)then 
     61            elseif(igas.eq.igas_C2H6)then
    6262               ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
    6363               mugaz_c = mugaz_c + 30.07*gfrac(igas)
    … …  
    9292         else
    9393            ! all values at 300 K from Engineering Toolbox
    94             if(igas.eq.igas_N2)then
     94            if(igas.eq.igas_CO2)then
    9595               !cpp_c   = cpp_c   + 0.744*gfrac(igas) ! @ ~210 K (better for
    96                !Mars conditions) 
     96               !Mars conditions)
    9797               cpp_c   = cpp_c   + 0.846*gfrac(igas)*44.01/mugaz_c
    9898            elseif(igas.eq.igas_N2)then

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

    r3184 r3353  
    22
    33!==================================================================
    4 !     
     4!
    55!     Purpose
    66!     -------
    7 !     Average the Rayleigh scattering in each band, weighting the 
     7!     Average the Rayleigh scattering in each band, weighting the
    88!     average by the blackbody function at temperature tstellar.
    9 !     Works for an arbitrary mix of gases (currently N2, N2 and 
     9!     Works for an arbitrary mix of gases (currently N2, N2 and
    1010!     H2 are possible).
    11 !     
     11!
    1212!     Authors
    13 !     ------- 
     13!     -------
    1414!     Robin Wordsworth (2010)
    1515!     Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995).
    … …  
    1919!     ---------
    2020!     setspv.F
    21 !     
     21!
    2222!     Calls
    2323!     -----
    2424!     none
    25 !     
     25!
    2626!==================================================================
    2727
    2828      use radinc_h, only: L_NSPECTV
    2929      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep
    30       use gases_h, only: ngasmx, vgas, gnom, gfrac, igas_N2, igas_H2, &
     30      use gases_h, only: ngasmx, vgas, gnom, gfrac, igas_CO2, igas_H2, &
    3131                           igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO
    3232      use comcstfi_mod, only: g, mugaz, pi
    … …  
    3737      integer N,Nfine,ifine,igas
    3838      parameter(Nfine=500.0)
    39       real*8 :: P0 = 9.423D+6   ! Rayleigh scattering reference pressure in pascals. P0 = 93*101325 Pa   
     39      real*8 :: P0 = 9.423D+6   ! Rayleigh scattering reference pressure in pascals. P0 = 93*101325 Pa
    4040
    4141      logical typeknown
    … …  
    7070         if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
    7171            print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
    72             'as its mixing ratio is less than 0.05.' 
     72            'as its mixing ratio is less than 0.05.'
    7373            ! ignore variable gas in Rayleigh calculation
    7474            ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
    7575            tauconsti(igas) = 0.0
    7676         else
    77             if(igas.eq.igas_N2) then
     77            if(igas.eq.igas_CO2) then
    7878               ! Hansen 1974 : equation (2.32)
    7979               tauconsti(igas) = (8.7/g)*1.527*scalep/P0
    … …  
    8282               tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
    8383            elseif(igas.eq.igas_H2O)then
    84                ! tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 
     84               ! tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0
    8585               tauconsti(igas) = 1.98017E-10*scalep/(g*mugaz*1E4)
    8686            elseif(igas.eq.igas_H2)then
    … …  
    103103               ! 16.04*1.6726*1E-27 is CH4 molecular mass
    104104               tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep
    105            
     105
    106106            elseif(igas.eq.igas_CO)then
    107107               ! see above for explanation
    … …  
    127127      endif
    128128
    129    
     129
    130130      do N=1,L_NSPECTV
    131          
     131
    132132         tausum = 0.0
    133133         tauwei = 0.0
    … …  
    147147                  tauvari(igas)   = 0.0
    148148               else
    149                   if(igas.eq.igas_N2)then
     149                  if(igas.eq.igas_CO2)then
    150150                     tauvari(igas) = (1.0+0.013/wl**2)/wl**4
    151151                  elseif(igas.eq.igas_N2)then
    … …  
    153153                  elseif(igas.eq.igas_H2O)then
    154154                     ! tauvari(igas) = 1.0/wl**4 ! to be improved...
    155                      tauvari(igas) = (5.7918E6/(2.38E2-1/wl**2)+1.679E5/(57.36E0-1/wl**2))**2/wl**4 
     155                     tauvari(igas) = (5.7918E6/(2.38E2-1/wl**2)+1.679E5/(57.36E0-1/wl**2))**2/wl**4
    156156                  elseif(igas.eq.igas_H2)then
    157                      tauvari(igas) = 1.0/wl**4 
     157                     tauvari(igas) = 1.0/wl**4
    158158
    159159                     ! below is exo_k formula : this tauvari is consistent with commented tauconsti for H2 above
    160160                     ! tauvari(igas) = (8.14E-49+1.28E-50/wl**2+1.61E-51/wl**4) / wl**4
    161161                  elseif(igas.eq.igas_He)then
    162                      tauvari(igas) = 1.0/wl**4 
     162                     tauvari(igas) = 1.0/wl**4
    163163                  elseif(igas.eq.igas_CH4)then
    164164                     tauvari(igas) = (4.6662E-4+4.02E-6/wl**2)**2 / wl**4
    … …  
    186186            tauweivar=tauweivar+df
    187187            tausumvar=tausumvar+tauvarvar*df
    188          
     188
    189189         enddo
    190190         TAURAY(N)=tausum/tauwei

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

    r3329 r3353  
    345345         CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,  &
    346346                                  pt,dtlw_hcn_c2h2)
    347         ! write(*,*) "Not supported yet"
    348         STOP
    349347      endif
    350348

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

    r3258 r3353  
    197197      logical,save :: oldplutocorrk
    198198!$OMP THREADPRIVATE(oldplutocorrk)
    199 
     199      logical,save :: oldplutosedim
     200!$OMP THREADPRIVATE(oldplutosedim)
    200201      logical,save :: global1d
    201202      real,save    :: szangle

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

    r3329 r3353  
    2525      character(len=300),save :: hazerad_file
    2626      character(len=300),save :: hazemmr_file
     27      character(len=300),save :: hazedens_file
     28
    2729
    2830      end module datafile_mod

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

    r3329 r3353  
    1313  use radcommon_h, only: ini_radcommon_h
    1414  use radii_mod, only: radfixed, Nmix_n2
    15   use datafile_mod, only: datadir, hazeprop_file, hazerad_file, hazemmr_file
     15  use datafile_mod, only: datadir,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file
    1616  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
    1717  use comgeomfi_h, only: totarea, totarea_planet
    … …  
    148148
    149149     if (is_master) write(*,*) "Haze mmr datafile"
    150      hazemmr_file="hazemmr.txt"  ! default file
     150     hazemmr_file="None"  ! default file
    151151     call getin_p("hazemmr_file",hazemmr_file)
    152152     if (is_master) write(*,*) trim(rname)//" hazemmr_file = ",trim(hazemmr_file)
    153 
     153     if (is_master) write(*,*) "Haze density datafile"
     154     hazedens_file="None"  ! default file
     155     call getin_p("hazedens_file",hazedens_file)
     156     if (is_master) write(*,*) trim(rname)//" hazedens_file = ",trim(hazedens_file)
    154157
    155158     if (is_master) write(*,*) trim(rname)//&
    … …  
    312315     call getin_p("UseTurbDiff",UseTurbDiff)
    313316     if (is_master) write(*,*) trim(rname)//": UseTurbDiff = ",UseTurbDiff
    314 
    315      if (is_master) write(*,*) trim(rname)//": use vdifc from old Pluto ?"
    316      oldplutovdifc=.true. ! default value
    317      call getin_p("oldplutovdifc",oldplutovdifc)
    318      if (is_master) write(*,*) trim(rname)//": oldplutovdifc = ",oldplutovdifc
    319 
    320      if (is_master) write(*,*) trim(rname)//&
    321        ": call pluto.old correlated-k radiative transfer ?"
    322      oldplutocorrk=.true. ! default value
    323      call getin_p("oldplutocorrk",oldplutocorrk)
    324      if (is_master) write(*,*) trim(rname)//": oldplutocorrk = ",oldplutocorrk
    325 
    326317
    327318     if (is_master) write(*,*) trim(rname)//": call convective adjustment ?"
    … …  
    636627      " thresh_non2 = ",thresh_non2
    637628
     629    ! use Pluto.old routines
     630     if (is_master) write(*,*) trim(rname)//": use vdifc from old Pluto ?"
     631     oldplutovdifc=.true. ! default value
     632     call getin_p("oldplutovdifc",oldplutovdifc)
     633     if (is_master) write(*,*) trim(rname)//": oldplutovdifc = ",oldplutovdifc
     634
     635     if (is_master) write(*,*) trim(rname)//&
     636       ": call pluto.old correlated-k radiative transfer ?"
     637     oldplutocorrk=.true. ! default value
     638     call getin_p("oldplutocorrk",oldplutocorrk)
     639     if (is_master) write(*,*) trim(rname)//": oldplutocorrk = ",oldplutocorrk
     640
     641     if (is_master) write(*,*) trim(rname)//&
     642       ": call pluto.old sedimentation ?"
     643     oldplutosedim=.true. ! default value
     644     call getin_p("oldplutosedim",oldplutosedim)
     645     if (is_master) write(*,*) trim(rname)//": oldplutosedim = ",oldplutosedim
     646
    638647!***************************************************************
    639648     !! Haze options
    … …  
    707716     if (is_master)write(*,*)trim(rname)//&
    708717     "nb_monomer = ",nb_monomer
     718
     719     if (is_master)write(*,*)trim(rname)//&
     720     "fixed gaseous CH4 mixing ratio for rad transfer ?"
     721    ch4fix=.false. ! default value
     722    call getin_p("ch4fix",ch4fix)
     723    if (is_master)write(*,*)trim(rname)//&
     724     " ch4fix = ",ch4fix
     725    if (is_master)write(*,*)trim(rname)//&
     726     "fixed gaseous CH4 mixing ratio for rad transfer ?"
     727    vmrch4fix=0.5 ! default value
     728    call getin_p("vmrch4fix",vmrch4fix)
     729    if (is_master)write(*,*)trim(rname)//&
     730     " vmrch4fix = ",vmrch4fix
     731    vmrch4_proffix=.false. ! default value
     732    call getin_p("vmrch4_proffix",vmrch4_proffix)
     733    if (is_master)write(*,*)trim(rname)//&
     734     " vmrch4_proffix = ",vmrch4_proffix
     735
     736    if (is_master)write(*,*)trim(rname)//&
     737     "call specific cooling for radiative transfer ?"
     738    cooling=.false.  ! default value
     739    call getin_p("cooling",cooling)
     740    if (is_master)write(*,*)trim(rname)//&
     741     " cooling = ",cooling
     742
     743    if (is_master)write(*,*)trim(rname)//&
     744     "NLTE correction for methane heating rates?"
     745    nlte=.false.  ! default value
     746    call getin_p("nlte",nlte)
     747    if (is_master)write(*,*)trim(rname)//&
     748     " nlte = ",nlte
     749    strobel=.false.  ! default value
     750    call getin_p("strobel",strobel)
     751    if (is_master)write(*,*)trim(rname)//&
     752     " strobel = ",strobel
    709753
    710754     if (is_master)write(*,*)trim(rname)//&

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

    r3329 r3353  
    5555                              n2cond,nearn2cond,noseason_day,conservn2, &
    5656                              convergeps,kbo,triton,paleo,paleoyears, &
    57                               carbox, methane, oldplutovdifc, oldplutocorrk, &
    58                               aerohaze,haze_proffix,source_haze, &
     57                              carbox, methane,&
     58                              oldplutovdifc,oldplutocorrk,oldplutosedim, &
     59                              aerohaze,haze_proffix,source_haze,&
    5960                              season, sedimentation,generic_condensation, &
    6061                              specOLR, &
    … …  
    979980                               cloudfrac,totcloudfrac,.false.,                   &
    980981                               firstcall,lastcall)
     982                  albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
    981983               else
    982984                call callcorrk(ngrid,nlayer,pq,nq,qsurf,  &
    … …  
    14801482            end if
    14811483
    1482             ! if (.not. water) then
    1483                ! Compute GCS (Generic Condensable Specie) cloud fraction. For now we can not have both water cloud fraction and GCS cloud fraction
    1484                ! Water is the priority
    1485                ! If you have set water and generic_condensation, then cloudfrac will be water cloudfrac
    1486                !
    1487                ! If you have set generic_condensation (and not water) and you have set several GCS
    1488                ! then cloudfrac will be only the cloudfrac of the last generic tracer
    1489                ! (Because it is rewritten every tracer in the loop)
    1490                !
    1491                ! Maybe one should create a cloudfrac_generic(ngrid,nlayer,nq) with 3 dimensions, the last one for tracers
    1492 
    14931484               ! Let's loop on tracers
    14941485               cloudfrac(:,:)=0.0
    … …  
    15071498         endif !generic_condensation
    15081499
    1509          !Generic Rain !AF24: removed
    1510 
    15111500         ! Sedimentation.
    15121501         if (sedimentation) then
    … …  
    15151504            zdqssed(1:ngrid,1:nq)  = 0.0
    15161505
    1517             ! if(watertest)then !AF24: removed
    1518 
    1519             call callsedim(ngrid,nlayer,ptimestep,           &
     1506            if (oldplutosedim)then
     1507               call callsedim_pluto(ngrid,nlayer,ptimestep,    &
     1508                          pplev,zzlev,pt,pdt,rice_ch4,rice_co, &
     1509                          pq,pdq,zdqsed,zdqssed,nq,pphi)
     1510
     1511            else
     1512               call callsedim(ngrid,nlayer,ptimestep,       &
    15201513                          pplev,zzlev,pt,pdt,pq,pdq,        &
    15211514                          zdqsed,zdqssed,nq)
    1522 
    1523             ! if(watertest)then !AF24: removed
     1515            endif
    15241516
    15251517            ! Whether it falls as rain or snow depends only on the surface temperature
    … …  
    15271519            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
    15281520
    1529 !!            call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
    1530 
    1531             ! ! Test water conservation !AF24: removed
    1532 
    15331521         end if ! end of 'sedimentation'
    1534 
    1535 !!         call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
    1536 !!         call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
    15371522
    15381523  ! ---------------
    … …  
    20962081         !call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
    20972082         call wstats(ngrid,"p","Pressure","Pa",3,pplay)
     2083         call wstats(ngrid,"emis","Emissivity","",2,emis)
    20982084         call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
    20992085         call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
    … …  
    21572143         endif
    21582144
     2145       call writediagfi(ngrid,"emis","Emissivity","",2,emis)
    21592146       call writediagfi(ngrid,"temp","temperature","K",3,zt)
    21602147       call writediagfi(ngrid,"teta","potential temperature","K",3,zh)

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

    r3329 r3353  
    106106         enddo
    107107         close(223)
    108          print*, 'TB22 READ HAZERAD'
     108         print*, 'Read hazerad from ',file_path
    109109       ENDIF
    110110
Note: See TracChangeset for help on using the changeset viewer.