Ignore:
Timestamp:
Apr 7, 2022, 4:33:05 PM (3 years ago)
Author:
cmathe
Message:

correction on rsedcloudco2/opacities/riceco2 computations, and more comments

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/co2cloud.F90

    r2616 r2660  
    6565!     3.2. No sub-grid cloud representation (CLFvaryingCO2 = False)
    6666!   4. Microphysics of CO2 cloud formation
    67 !     4.1. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
    68 !     4.2. Effective tracers quantities in the cloud
    69 !     4.3. Gravitational sedimentation
    70 !       4.3.a. Compute cloud density
    71 !       4.3.b. Save actualized tracer values to compute sedimentation tendancies
    72 !       4.3.c. Sedimentation of co2 ice
    73 !       4.3.d. Sedimentation for other tracers
    74 !       4.3.e. Compute tendencies due to the sedimation process
    75 !     4.4. Main call to the cloud scheme
    76 !     4.5. Updating tendencies after cloud scheme
     67!     4.1. Effective tracers quantities in the cloud
     68!     4.2. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
     69!     4.3. Main call to the cloud scheme
     70!     4.4. Updating tendencies after cloud scheme
     71!     4.5. Gravitational sedimentation
     72!       4.5.a. Compute cloud density
     73!       4.5.b. Save actualized tracer values to compute sedimentation tendancies
     74!       4.5.c. Sedimentation of co2 ice
     75!       4.5.d. Sedimentation for other tracers
     76!       4.5.e. Compute tendencies due to the sedimation process
    7777!   5. Compute final tendencies after time loop
    7878!   6. Update clouds physical values in the cloud (for output)
    … …  
    8080!     6.2. Update rice and rdust
    8181!   7. Correction if a lot of subliming CO2 fills the 1st layer
    82 !   8. Compute water cloud sedimentation radius
    83 !   9. CO2 saturated quantities
    84 !     9.1 Compute CO2 saturation in layers
    85 !     9.2 Compute CO2 saturated quantities in layers
    86 !   10. Everything modified by CO2 microphysics must be wrt co2cloudfrac
    87 !   11. Compute opacity at 1 micron
    88 !   12. Write outputs in diagfi.nc
     82!   8. CO2 saturated quantities
     83!     8.1 Compute CO2 saturation in layers
     84!     8.2 Compute CO2 saturated quantities in layers
     85!   9. Everything modified by CO2 microphysics must be wrt co2cloudfrac
     86!   10. Compute opacity at 1 micron
     87!   11. Write outputs in diagfi.nc
    8988!======================================================================================================================!
    9089  subroutine co2cloud(ngrid, nlay, ptimestep, pplev, pplay, pdpsrf, pzlay, pt, pdt, pq, pdq, pdqcloudco2, pdtcloudco2, &
    … …  
    106105
    107106  use datafile_mod,    only: datadir
     107  use density_co2_ice_mod, only: density_co2_ice
    108108
    109109  use improvedCO2clouds_mod, only: improvedCO2clouds
    … …  
    173173
    174174  real, parameter :: &
    175      mincloud = 0.1,  &! Minimum cloud fraction
    176      beta = 0.85,     &! correction for the shape of the particles (see Murphy et al. JGR 1990 vol.95):
    177                        !   beta = 1    for spheres
    178                        !   beta = 0.85 for irregular particles
    179                        !   beta = 0.5  for disk shaped particles
    180      threshold = 1e-30 ! limit value
     175     mincloud = 0.1,     &! Minimum cloud fraction
     176     beta = 0.85,        &! correction for the shape of the particles (see Murphy et al. JGR 1990 vol.95):
     177                          !   beta = 1    for spheres
     178                          !   beta = 0.85 for irregular particles
     179                          !   beta = 0.5  for disk shaped particles
     180     threshold = 1e-30, & ! limit value considering as zero
     181     threshold_2 = 1e-13  ! limit value considering the value is physical (below this value => computer noise for dble)
    181182
    182183  double precision, parameter :: &
    … …  
    550551    ! TODO: Totalcloud frac of the column missing here
    551552!----------------------------------------------------------------------------------------------------------------------!
    552 ! 3.2. No sub-grid cloud representation (CLFvarying = False)
     553! 3.2. No sub-grid cloud representation (CLFvaryingCO2 = False)
    553554!----------------------------------------------------------------------------------------------------------------------!
    554555  else
    … …  
    565566  pteff(:,:) = pt(:,:)
    566567!----------------------------------------------------------------------------------------------------------------------!
    567 ! 4.2. Effective tracers quantities in the cloud
     568! 4.1. Effective tracers quantities in the cloud
    568569!----------------------------------------------------------------------------------------------------------------------!
    569570  if (CLFvaryingCO2) then
    … …  
    577578  do microstep = 1, imicroco2
    578579!----------------------------------------------------------------------------------------------------------------------!
    579 ! 4.1. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
     580! 4.2. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
    580581!----------------------------------------------------------------------------------------------------------------------!
    581582   do l = 1, nlay
    … …  
    623624   end do ! nlay
    624625!----------------------------------------------------------------------------------------------------------------------!
    625 ! 4.4. Main call to the cloud scheme
     626! 4.3. Main call to the cloud scheme
    626627!----------------------------------------------------------------------------------------------------------------------!
    627628   call improvedco2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, &
    … …  
    695696   end do
    696697!----------------------------------------------------------------------------------------------------------------------!
    697 ! 4.5. Updating tendencies after cloud scheme
     698! 4.4. Updating tendencies after cloud scheme
    698699!----------------------------------------------------------------------------------------------------------------------!
    699700    do l = 1, nlay
    … …  
    741742    end do ! nlay
    742743!----------------------------------------------------------------------------------------------------------------------!
    743 ! 4.3. Gravitational sedimentation
     744! 4.5. Gravitational sedimentation
    744745!----------------------------------------------------------------------------------------------------------------------!
    745746    if (sedimentation) then
    746747!----------------------------------------------------------------------------------------------------------------------!
    747 ! 4.3.a. Compute cloud density
     748! 4.5.a. Compute cloud density
    748749!----------------------------------------------------------------------------------------------------------------------!
    749750      do l = 1, nlay
    … …  
    768769            Nccnco2 = Nccnco2 - Nccnco2_h2o
    769770            Qccnco2 = Qccnco2 - Qccnco2_h2o
     771            if (Nccnco2 < 0.) then
     772              Nccnco2 = threshold
     773              Qccnco2 = threshold
     774            end if
    770775          end if
    771776
    … …  
    781786      end do ! nlay
    782787!----------------------------------------------------------------------------------------------------------------------!
    783 ! 4.3.b. Save actualized tracer values to compute sedimentation tendancies
     788! 4.5.b. Save actualized tracer values to compute sedimentation tendancies
    784789!----------------------------------------------------------------------------------------------------------------------!
    785790      zqsed0(:,:,igcm_co2_ice) = zqsed(:,:,igcm_co2_ice)
    … …  
    798803      end if
    799804!----------------------------------------------------------------------------------------------------------------------!
    800 ! 4.3.c. Sedimentation of co2 ice
     805! 4.5.c. Sedimentation of co2 ice
    801806!----------------------------------------------------------------------------------------------------------------------!
    802807      do ig = 1, ngrid
    803808        do l = 1, nlay
    804           rsedcloudco2(ig,l) = max( riceco2(ig,l)*(1.+sigma_iceco2)*(1.+sigma_iceco2)*(1.+sigma_iceco2), &
     809          rsedcloudco2(ig,l) = max( riceco2(ig,l)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed), &
    805810                                    rdust(ig,l) )
    806811        end do
    … …  
    815820      end do
    816821!----------------------------------------------------------------------------------------------------------------------!
    817 ! 4.3.d. Sedimentation for other tracers
     822! 4.5.d. Sedimentation for other tracers
    818823!----------------------------------------------------------------------------------------------------------------------!
    819824      wq(:,:) = 0.
    … …  
    881886      end if
    882887!----------------------------------------------------------------------------------------------------------------------!
    883 ! 4.3.e. Compute tendencies due to the sedimation process
     888! 4.5.e. Compute tendencies due to the sedimation process
    884889!----------------------------------------------------------------------------------------------------------------------!
    885890      do l = 1, nlay
    … …  
    961966    do ig = 1, ngrid
    962967      pdtcloudco2(ig,l) = ( sum_subpdt(ig,l)/real(imicroco2) ) - pdt(ig,l)
     968      if (pdtcloudco2(ig,l) /= pdtcloudco2(ig,l)) then
     969        write(*,*)ig,l,pdtcloudco2(ig,l), pdt(ig,l), sum_subpdt(ig,l)
     970        call abort_physic('co2clouds', 'ptdcloudco2 is NaN', 1)
     971      end if
    963972    end do
    964973  end do
    … …  
    967976  do l = 1, nlay
    968977    do ig = 1, ngrid
    969       pdqcloudco2(ig,l,igcm_co2) = 0. ! here is the trick, this tendency is computed in co2condens_mod4micro
     978      pdqcloudco2(ig,l,igcm_co2) = 0. ! here is the trick, this tendency is computed in co2condens
    970979
    971980      pdqcloudco2(ig,l,igcm_co2_ice) = ( sum_subpdq(ig,l,igcm_co2_ice) / real(imicroco2) ) - pdq(ig,l,igcm_co2_ice)
    … …  
    10201029    do ig = 1, ngrid
    10211030      Niceco2 = pqeff(ig,l,igcm_co2_ice) + ( pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice) ) * ptimestep
    1022       Niceco2 = max (Niceco2, threshold)
     1031      Niceco2 = max(Niceco2, threshold)
    10231032
    10241033      ! meteoritic particles are considered like dust, => rho_dust
    … …  
    10471056      Nccnco2 = Nccnco2 - Nccnco2_h2o
    10481057      Qccnco2 = Qccnco2 - Qccnco2_h2o
     1058        if (Nccnco2 <= 0) then
     1059          Nccnco2 = threshold
     1060          Qccnco2 = threshold
     1061        end if
    10491062      end if
     1063
    10501064      ! Compute particle size
    10511065      call updaterice_microco2(dble(Niceco2), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), dble(Nccnco2_h2o), myT, &
    … …  
    10531067
    10541068     ! Compute opacities
    1055       if ( (Niceco2 <= threshold .or. (Nccnco2)*tauscaling(ig) <= 1.) ) then ! Nccnco2 inclut Nccnco2_h2o
     1069      if ( (Niceco2 >= threshold_2 .and. (Nccnco2+Nccnco2_h2o)*tauscaling(ig) >= 1.) ) then
     1070        n_derf = derf( (rb_cldco2(1)-log(riceco2(ig,l))) *dev2)
     1071        Qext1bins2(ig,l) = 0.
     1072        do i = 1, nbinco2_cld
     1073          n_aer(i) = -0.5 * (Nccnco2+Nccnco2_h2o)*tauscaling(ig) * n_derf
     1074
     1075          n_derf = derf((rb_cldco2(i+1)-log(riceco2(ig,l))) *dev2)
     1076          n_aer(i) = n_aer(i) + (0.5 * (Nccnco2+Nccnco2_h2o)*tauscaling(ig) * n_derf)
     1077
     1078          Qext1bins2(ig,l) = Qext1bins2(ig,l) + Qext1bins(i) * n_aer(i)
     1079        end do
     1080      else
    10561081        riceco2(ig,l) = 0.
    10571082        Qext1bins2(ig,l) = 0.
    1058       else
    1059         n_derf = derf( (rb_cldco2(1)-log(riceco2(ig,l))) *dev2)
    1060         Qext1bins2(ig,l) = 0.
    1061 
    1062         do i = 1, nbinco2_cld
    1063           n_aer(i) = -0.5 * (Nccnco2)*tauscaling(ig) * n_derf
    1064 
    1065           n_derf = derf((rb_cldco2(i+1)-log(riceco2(ig,l))) *dev2)
    1066           n_aer(i) = n_aer(i) + (0.5 * (Nccnco2)*tauscaling(ig) * n_derf)
    1067 
    1068           Qext1bins2(ig,l) = Qext1bins2(ig,l) + Qext1bins(i) * n_aer(i)
    1069         end do
    10701083      end if
    10711084!----------------------------------------------------------------------------------------------------------------------!
    … …  
    11031116  end do
    11041117!----------------------------------------------------------------------------------------------------------------------!
    1105 ! 9. CO2 saturated quantities
    1106 !----------------------------------------------------------------------------------------------------------------------!
    1107 ! 9.1 Compute CO2 saturation in layers
     1118! 8. CO2 saturated quantities
     1119!----------------------------------------------------------------------------------------------------------------------!
     1120! 8.1 Compute CO2 saturation in layers
    11081121!----------------------------------------------------------------------------------------------------------------------!
    11091122  call co2sat(ngrid*nlay, pteff+(pdt+pdtcloudco2)*ptimestep, zqsatco2)
    11101123!----------------------------------------------------------------------------------------------------------------------!
    1111 ! 9.2 Compute CO2 saturated quantities in layers
     1124! 8.2 Compute CO2 saturated quantities in layers
    11121125!----------------------------------------------------------------------------------------------------------------------!
    11131126  do l = 1, nlay
    … …  
    11181131  end do
    11191132!----------------------------------------------------------------------------------------------------------------------!
    1120 ! 10. Everything modified by CO2 microphysics must be wrt co2cloudfrac
     1133! 9. Everything modified by CO2 microphysics must be wrt co2cloudfrac
    11211134!----------------------------------------------------------------------------------------------------------------------!
    11221135  if (CLFvaryingCO2) then
    … …  
    11611174  end if ! if CLFvaryingCO2 is true
    11621175!----------------------------------------------------------------------------------------------------------------------!
    1163 ! 11. Compute opacity at 1 micron: Opacity in mesh ig is the sum over l of Qext1bins2. Is this true ?
     1176! 10. Compute opacity at 1 micron: Opacity in mesh ig is the sum over l of Qext1bins2. Is this true ?
    11641177!----------------------------------------------------------------------------------------------------------------------!
    11651178  tau1mic(:)=0.
    … …  
    11701183  end do
    11711184!----------------------------------------------------------------------------------------------------------------------!
    1172 ! 12. Write outputs in diagfi.nc
     1185! 11. Write outputs in diagfi.nc
    11731186!----------------------------------------------------------------------------------------------------------------------!
    11741187  call WRITEDIAGFI(ngrid, "satuco2", "vap in satu", " ", 3, satuco2)

  • trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F

    r2617 r2660  
    335335     &                       rdust,zcondicea,zfallice,zdq_scav,pdqsc)
    336336      endif
    337      
     337             call WRITEdiagfi(ngrid,"zfallice",
     338     &         "",
     339     &         " ",2,zfallice(ngrid,1))
    338340      ELSE ! if co2 clouds
    339341        condens_layer(:,:) = 0.
    … …  
    350352          piceco2(ig) = piceco2(ig) + zdqssed_co2(ig)*ptimestep
    351353        ENDDO
     354       call WRITEdiagfi(ngrid,"zfallice",
     355     &         "",
     356     &         " ",2,zdqssed_co2) ! otherwise we have not 1 day(1proc) = 1 day (24procs) test
    352357      ENDIF ! end of if co2clouds
    353358
    … …  
    359364     &         "",
    360365     &         " ",3,condens_layer)
    361      
    362        call WRITEdiagfi(ngrid,"zfallice",
    363      &         "",
    364      &         " ",2,zfallice(ngrid,1))
    365      
    366366
    367367c     *************************
    … …  
    441441c           """"""""""""""""""""""""""""""""""""""""""""""""""""
    442442c           (including what has just condensed in the atmosphere)
     443            if (co2clouds) then
     444            IF((piceco2(ig)/ptimestep).LE.
     445     &          -zcondices(ig))THEN
     446               zcondices(ig) = -piceco2(ig)/ptimestep
     447               pdtsrfc(ig)=(latcond/pcapcal(ig)) *
     448     &         (zcondices(ig)+zfallheat)
     449            END IF
     450            else
    443451            IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
    444452     &          -zcondices(ig))THEN
    445453               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
    446                pdtsrfc(ig)=(latcond/pcapcal(ig))*
     454               pdtsrfc(ig)=(latcond/pcapcal(ig)) *
    447455     &         (zcondices(ig)+zfallheat)
    448456            END IF
     457            end if
    449458
    450459c           Changing CO2 ice amount and pressure :

  • trunk/LMDZ.MARS/libf/phymars/improvedco2clouds_mod.F90

    r2616 r2660  
    133133
    134134  real, parameter :: &
    135      threshold = 1e-30 ! limit value
    136 
    137   character(len=23), parameter:: &
     135     threshold = 1e-30, & ! limit value considering as zero
     136     threshold_2 = 1e-13  ! limit value considering the value is physical (below this value => computer noise for dble)
     137
     138  character(len=50), parameter:: &
    138139     file_meteoritic_flux = 'Meteo_flux_Plane.dat'
    139140!----------------------------------------------------------------------------------------------------------------------!
    … …  
    194195    pco2,                              &! Co2 vapor partial pressure (Pa)
    195196    satu,                              &! Co2 vapor saturation ratio over ice
    196     Mo,                                &! mass of aerosol particles
    197     No,                                &! number of aerosol particles
     197    Mo,                                &! mass of dust particles
     198    No,                                &! number of dust particles
    198199    Rn,                                &! logarithm of rdust/rice
    199200    Rm,                                &! Rn * variance of ice and CCN distribution
    200201    n_derf,                            &! derf( (rb_cldco2(1)+Rn) *dev3)
    201202    m_derf,                            &! derf( (rb_cldco2(1)+Rm) *dev2)
    202     n_aer(nbinco2_cld),                &! Radius used by the microphysical scheme (m)
    203     m_aer(nbinco2_cld),                &! number concentration V-1 of particle/each size bin
    204     n_aer_meteor(nbinco2_cld),                &! Radius used by the microphysical scheme (m)
    205     m_aer_meteor(nbinco2_cld),                &! number concentration V-1 of particle/each size bin
    206     n_aer_h2oice(nbinco2_cld),         &! mass mixing ratio of particle/each size bin
    207     m_aer_h2oice(nbinco2_cld),         &! Same - for CO2 nucleation
     203    n_aer(nbinco2_cld),                &! Radius used by the microphysical scheme (m) for dust
     204    m_aer(nbinco2_cld),                &! number concentration V-1 of particle/each size bin for dust
     205    n_aer_meteor(nbinco2_cld),         &! Radius used by the microphysical scheme (m) for meteor
     206    m_aer_meteor(nbinco2_cld),         &! number concentration V-1 of particle/each size bin for meteor
     207    n_aer_h2oice(nbinco2_cld),         &! mass mixing ratio of particle/each size bin for h2o ice
     208    m_aer_h2oice(nbinco2_cld),         &! Same - for CO2 nucleation for h2o ice
    208209    Ic_rice,                           &! Mass transfer rate CO2 ice crystal
    209210    ratioh2o_ccn,                      &! 1./(zq(ig,l,igcm_h2o_ice)  + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
    … …  
    211212    dN,                                &! number of particle of dust used as ccn
    212213    dM,                                &! mass of dN
    213     dN_meteor,                                &! number of particle of dust used as ccn
    214     dM_meteor,                                &! mass of dN
     214    dN_meteor,                         &! number of particle of meteoritic particles used as ccn
     215    dM_meteor,                         &! mass of dN_meteor
    215216    dNh2o,                             &! number of particle of h2o ice used as ccn
    216217    dMh2o,                             &! mass of dNh2o
    217218    dNN,                               &! min(dN,zq(ig,l,igcm_dust_number))
    218219    dMM,                               &! min(dM,zq(ig,l,igcm_dust_mass))
    219     dNN_meteor,                                &! number of particle of dust used as ccn
    220     dMM_meteor,                                &! mass of dN
    221220    dNNh2o,                            &! min(dNNh2o,zq(ig,l,igcm_ccn_number))
    222221    dMh2o_ice,                         &! min(dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn, zq(ig,l,igcm_h2o_ice))
    223222    dMh2o_ccn,                         &! min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
    224223    rate(nbinco2_cld),                 &! nucleation rate
    225     rate_meteor(nbinco2_cld),                 &! nucleation rate
     224    rate_meteor(nbinco2_cld),          &! nucleation rate for meteor
    226225    rateh2o(nbinco2_cld),              &! nucleation rate for h2o
    227226    rho_ice_co2T,                      &! density of co2 ice Temperature-dependent
    … …  
    230229    Proba,                             &! 1.0 - exp(-1.*microtimestep*rate(i))
    231230    Probah2o,                          &! 1.0 - exp(-1.*microtimestep*rateh2o(i))
    232     Proba_meteor,                          &! 1.0 - exp(-1.*microtimestep*rateh2o(i))
     231    Proba_meteor,                      &! 1.0 - exp(-1.*microtimestep*rate_meteor(i))
    233232    diff_pressure(nlev_meteor),        &! abs(pression_meteor(:)-pplev(ig,l))
    234     meteor_ccn(ngrid,nlay,nbinco2_cld)  !
     233    meteor_ccn(ngrid,nlay,nbinco2_cld)  ! input flux of meteoritc particles
    235234
    236235  logical :: &
    … …  
    253252    if (meteo_flux) then
    254253      ! Check if file exists
    255       inquire(file=trim(datadir)//'/'//file_meteoritic_flux, exist=file_ok)
     254      inquire(file=trim(datadir)//'/'//trim(file_meteoritic_flux), exist=file_ok)
    256255      if (.not. file_ok) then
    257         call abort_physic("CO2clouds", 'file '//file_meteoritic_flux//' should be in'//trim(datadir), 1)
     256        call abort_physic("CO2clouds", 'file '//trim(file_meteoritic_flux)//' should be in'//trim(datadir), 1)
    258257      end if
     258      write(*,*)'Meteoritic flux file used: ', trim(file_meteoritic_flux)
    259259
    260260      ! open file
    261       open(unit=uMeteor,file=trim(datadir)//'/'//file_meteoritic_flux, FORM='formatted')
     261      open(unit=uMeteor,file=trim(datadir)//'/'//trim(file_meteoritic_flux), FORM='formatted')
    262262
    263263      !skip 1 line
    … …  
    509509          Nccnco2 = Nccnco2 - Nccnco2_h2o
    510510          Qccnco2 = Qccnco2 - Qccnco2_h2o
     511          if (Nccnco2 <= 0.) then
     512            Nccnco2 = threshold
     513            Qccnco2 = threshold
     514          end if
    511515        end if
    512516
    … …  
    550554!----------------------------------------------------------------------------------------------------------------------!
    551555      ! On sublime tout
    552       if ((zq(ig,l,igcm_co2_ice) <= threshold).or.(zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) < 1.)) then
     556      if ((zq(ig,l,igcm_co2_ice) < threshold_2).or.(zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) < 1.)) then
    553557          zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccnco2_mass)
    554558          zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccnco2_number)

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r2646 r2660  
    509509      REAL co2totA
    510510      REAL co2totB
     511      REAL co2conservation
    511512
    512513c entrainment by mountain top dust flows above sub-grid scale topography
    … …  
    20312032         zdvc(:,:) = 0.
    20322033         zdqc(:,:,:) = 0.
     2034         zdqsc(:,:) = 0.
    20332035         CALL co2condens(ngrid,nlayer,nq,ptimestep,
    20342036     $              capcal,zplay,zplev,tsurf,pt,
    … …  
    39763978        enddo
    39773979      endif ! of if (igcm_co2_ice.ne.0)
    3978         call WRITEDIAGFI(ngrid,'co2conservation',
     3980      co2conservation = (co2totA-co2totB)/co2totA
     3981        call WRITEDIAGFI(ngrid, 'co2conservation',
    39793982     &                     'Total CO2 mass conservation in physic',
    3980      &                     '%',0,[(co2totA-co2totB)/co2totA])
     3983     &                     ' ', 0, co2conservation)
    39813984! XIOS outputs
    39823985#ifdef CPP_XIOS     

  • trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F

    r2628 r2660  
    9090c     converted to part kg-1 using a typical atmospheric density.
    9191
    92       REAL, PARAMETER :: ccn0 = 1.3E8
     92      REAL, PARAMETER :: ccn0 = 1.3E8, threshold = 1e-30 ! limit value
    9393     
    9494c     For microphysics only:     
    … …  
    9999
    100100      LOGICAL,SAVE :: firstcall=.true.
    101       REAL Nccnco2, Qccnco2
     101      REAL Nccnco2, Qccnco2, Niceco2
    102102      REAL Nccnco2_h2o, Qccnco2_h2o
    103103      REAL CBRT
    … …  
    198198          DO l=1,nlayer
    199199            DO ig=1,ngrid
    200               Nccnco2 = pq(ig,l,igcm_ccnco2_number)
    201               Qccnco2 = pq(ig,l,igcm_ccnco2_mass)
     200              Niceco2 = max(pq(ig,l,igcm_co2_ice), threshold)
     201              Nccnco2 = max(pq(ig,l,igcm_ccnco2_number),
     202     &                              threshold)
     203              Qccnco2 = max(pq(ig,l,igcm_ccnco2_mass),
     204     &                              threshold)
    202205              Nccnco2_h2o = 0.
    203206              Qccnco2_h2o = 0.
    204207              if (co2useh2o) then
    205                 Nccnco2_h2o = pq(ig,l,igcm_ccnco2_h2o_number)
    206                 Qccnco2_h2o = pq(ig,l,igcm_ccnco2_h2o_mass_ice)
    207      &                        + pq(ig,l,igcm_ccnco2_h2o_mass_ccn)
     208                Nccnco2_h2o = max(pq(ig,l,igcm_ccnco2_h2o_number),
     209     &                              threshold)
     210                Qccnco2_h2o = max(pq(ig,l,igcm_ccnco2_h2o_mass_ice)
     211     &                        + pq(ig,l,igcm_ccnco2_h2o_mass_ccn),
     212     &                              threshold)
     213
    208214                Nccnco2 = Nccnco2 - Nccnco2_h2o
    209215                Qccnco2 = Qccnco2 - Qccnco2_h2o
     216                if (Nccnco2 <= 0) then
     217                  Nccnco2 = threshold
     218                  Qccnco2 = threshold
     219                end if
    210220              end if
    211 
    212               call updaterice_microco2(dble(pq(ig,l,igcm_co2_ice)),
     221              call updaterice_microco2(dble(Niceco2),
    213222     &                              dble(Qccnco2), dble(Nccnco2),
    214223     &                              dble(Qccnco2_h2o),
Note: See TracChangeset for help on using the changeset viewer.