Index: trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F90 (revision 4183) +++ trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F90 (revision 4186) @@ -9,18 +9,18 @@ SUBROUTINE improvedclouds(ngrid,nlay,ptimestep,pplay,pt,pdt,pq,pdq,nq,tauscaling,imicro,zt,zq) -use updaterad, only: updaterice_micro, updaterccn -use watersat_mod, only: watersat -use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap, igcm_h2o_ice, igcm_dust_mass, & - igcm_dust_number, igcm_ccn_mass, igcm_ccn_number, igcm_hdo_vap, & - igcm_hdo_ice,qperemin -use conc_mod, only: mmean -use comcstfi_h, only: pi, cpp -use microphys_h, only: nbin_cld, rad_cld, mteta, kbz, nav, rgp -use microphys_h, only: mco2, vo1, mh2o, mhdo, molco2, molhdo, To -use nuclea_mod, only: nuclea -use sig_h2o_mod, only: sig_h2o -use growthrate_mod, only: growthrate +use updaterad, only: updaterice_micro, updaterccn +use watersat_mod, only: watersat +use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap, igcm_h2o_ice, igcm_dust_mass, & + igcm_dust_number, igcm_ccn_mass, igcm_ccn_number, igcm_hdo_vap, & + igcm_hdo_ice, qperemin +use conc_mod, only: mmean +use comcstfi_h, only: pi, cpp +use microphys_h, only: nbin_cld, rad_cld, mteta, kbz, nav, rgp +use microphys_h, only: mco2, vo1, mh2o, mhdo, molco2, molhdo, To +use nuclea_mod, only: nuclea +use sig_h2o_mod, only: sig_h2o +use growthrate_mod, only: growthrate use write_output_mod, only: write_output -use callkeys_mod, only: activice, scavenging, cloud_adapt_ts, hdo, hdofrac +use callkeys_mod, only: activice, scavenging, cloud_adapt_ts, hdo, hdofrac implicit none @@ -50,78 +50,78 @@ !------------------------------------------------------------------ ! Inputs/outputs: -INTEGER, INTENT(IN) :: ngrid,nlay -INTEGER, INTENT(IN) :: nq ! nombre de traceurs -REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s) -REAL, dimension(ngrid,nlay), INTENT(IN) :: pplay ! pression au milieu des couches (Pa) -REAL, dimension(ngrid,nlay), INTENT(IN) :: pt ! temperature at the middle of the layers (K) -REAL, dimension(ngrid,nlay), INTENT(IN) :: pdt ! temperature tendency (K/s) -REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pq ! tracer (kg/kg) -REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pdq ! tracer tendency (kg/kg/s) -REAL, dimension(ngrid), INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust -INTEGER, INTENT(IN) :: imicro ! nb. microphy calls(retrocompatibility) - -REAL, dimension(ngrid,nlay,nq), INTENT(OUT) :: zq ! tracers post microphy (kg/kg) -REAL, dimension(ngrid,nlay), INTENT(OUT) :: zt ! temperature post microphy (K) +INTEGER, INTENT(IN) :: ngrid,nlay +INTEGER, INTENT(IN) :: nq ! nombre de traceurs +REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s) +REAL, dimension(ngrid,nlay), INTENT(IN) :: pplay ! pression au milieu des couches (Pa) +REAL, dimension(ngrid,nlay), INTENT(IN) :: pt ! temperature at the middle of the layers (K) +REAL, dimension(ngrid,nlay), INTENT(IN) :: pdt ! temperature tendency (K/s) +REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pq ! tracer (kg/kg) +REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pdq ! tracer tendency (kg/kg/s) +REAL, dimension(ngrid), INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust +INTEGER, INTENT(IN) :: imicro ! nb. microphy calls(retrocompatibility) + +REAL, dimension(ngrid,nlay,nq), INTENT(OUT) :: zq ! tracers post microphy (kg/kg) +REAL, dimension(ngrid,nlay), INTENT(OUT) :: zt ! temperature post microphy (K) !------------------------------------------------------------------ ! Local variables: -LOGICAL, SAVE :: firstcall = .true. +LOGICAL, SAVE :: firstcall = .true. !$OMP THREADPRIVATE(firstcall) -REAL*8 :: derf ! Error function +REAL*8 :: derf ! Error function !external derf -INTEGER :: ig,l,i -REAL, dimension(ngrid,nlay) :: zqsat ! saturation -REAL :: lw ! Latent heat of sublimation (J.kg-1) -REAL :: cste -REAL :: dMice ! mass of condensed ice -REAL :: dMicetot ! JN -REAL :: dMice_hdo ! mass of condensed HDO ice -REAL, dimension(ngrid,nlay) :: alpha ! HDO equilibrium fractionation coefficient (Saturation=1) -REAL, dimension(ngrid,nlay) :: alpha_c ! HDO real fractionation coefficient -REAL*8 :: ph2o ! Water vapor partial pressure (Pa) -REAL*8 :: satu ! Water vapor saturation ratio over ice -REAL*8 :: Mo,No, Rn, Rm, dev2, n_derf, m_derf -REAL*8, dimension(nbin_cld) :: n_aer ! number conc. of particle/each size bin -REAL*8, dimension(nbin_cld) :: m_aer ! mass mixing ratio of particle/each size bin +INTEGER :: ig, l, i +REAL, dimension(ngrid,nlay) :: zqsat ! Saturation +REAL :: lw ! Latent heat of sublimation (J.kg-1) +REAL :: cste +REAL :: dMice ! mass of condensed ice +REAL :: dMicetot +REAL :: dMice_hdo ! mass of condensed HDO ice +REAL, dimension(ngrid,nlay) :: alpha ! HDO equilibrium fractionation coefficient (Saturation=1) +REAL, dimension(ngrid,nlay) :: alpha_c ! HDO real fractionation coefficient +REAL*8 :: ph2o ! Water vapor partial pressure (Pa) +REAL*8 :: satu ! Water vapor saturation ratio over ice +REAL*8 :: Mo, No, Rn, Rm, dev2, n_derf, m_derf +REAL*8, dimension(nbin_cld) :: n_aer ! number conc. of particle/each size bin +REAL*8, dimension(nbin_cld) :: m_aer ! mass mixing ratio of particle/each size bin REAL :: dN, dM, seq -REAL, dimension(nbin_cld) :: rate ! nucleation rate -REAL, dimension(ngrid,nlay) :: rice ! Ice mass mean radius (m) (r_c in montmessin_2004) -REAL, dimension(ngrid,nlay) :: rhocloud ! Cloud density (kg.m-3) -REAL, dimension(ngrid,nlay) :: rdust ! Dust geometric mean radius (m) -REAL :: res ! Resistance growth -REAL :: Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient +REAL, dimension(nbin_cld) :: rate ! Nucleation rate +REAL, dimension(ngrid,nlay) :: rice ! Ice mass mean radius (m) (r_c in montmessin_2004) +REAL, dimension(ngrid,nlay) :: rhocloud ! Cloud density (kg.m-3) +REAL, dimension(ngrid,nlay) :: rdust ! Dust geometric mean radius (m) +REAL :: res ! Resistance growth +REAL :: Dv, Dv_hdo ! Water/HDO vapor diffusion coefficient ! Parameters of the size discretization used by the microphysical scheme -DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m) -DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m) -DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6 ! Minimum boundary radius (m) -DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m) -DOUBLE PRECISION :: vrat_cld ! Volume ratio -DOUBLE PRECISION, dimension(nbin_cld+1), SAVE :: rb_cld ! boundary values of each rad_cld bin (m) +DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m) +DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m) +DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6 ! Minimum boundary radius (m) +DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m) +DOUBLE PRECISION :: vrat_cld ! Volume ratio +DOUBLE PRECISION, dimension(nbin_cld + 1), SAVE :: rb_cld ! Boundary values of each rad_cld bin (m) !$OMP THREADPRIVATE(rb_cld) -DOUBLE PRECISION, dimension(nbin_cld) :: dr_cld ! width of each rad_cld bin (m) -DOUBLE PRECISION, dimension(nbin_cld) :: vol_cld ! particle volume for each bin (m3) -REAL, SAVE :: sigma_ice ! Variance of the ice and CCN distributions +DOUBLE PRECISION, dimension(nbin_cld) :: dr_cld ! Width of each rad_cld bin (m) +DOUBLE PRECISION, dimension(nbin_cld) :: vol_cld ! Particle volume for each bin (m3) +REAL, SAVE :: sigma_ice ! Variance of the ice and CCN distributions !$OMP THREADPRIVATE(sigma_ice) -!---------------------------------- -! JN : used in subtimestep -REAL :: microtimestep! subdivision of ptimestep (s) -REAL, dimension(ngrid,nlay) :: subpdtcloud ! Temperature variation due to latent heat -REAL, dimension(ngrid,nlay,nq) :: zq0 ! local initial value of tracers -INTEGER, dimension(ngrid,nlay) :: zimicro ! Subdivision of ptimestep -INTEGER, dimension(ngrid,nlay) :: count_micro ! Number of microphys calls -REAL, dimension(ngrid,nlay) :: zpotcond ! maximal condensable water (previous two) -REAL, dimension(1) :: zqsat_tmp ! maximal condensable water (previous two) -REAL :: spenttime ! time spent in while loop -REAL :: zdq ! used to compute adaptative timestep -LOGICAL :: ending_ts ! Condition to end while loop - -!---------------------------------- +!---------------------------------- +! JN: used in subtimestep +REAL :: microtimestep ! Subdivision of ptimestep (s) +REAL, dimension(ngrid,nlay) :: subpdtcloud ! Temperature variation due to latent heat +REAL, dimension(ngrid,nlay,nq) :: zq0 ! Local initial value of tracers +INTEGER, dimension(ngrid,nlay) :: zimicro ! Subdivision of ptimestep +INTEGER, dimension(ngrid,nlay) :: count_micro ! Number of microphys calls +REAL, dimension(ngrid,nlay) :: zpotcond ! Maximal condensable water (previous two) +REAL, dimension(1) :: zqsat_tmp ! Maximal condensable water (previous two) +REAL :: spenttime ! Time spent in while loop +REAL :: zdq ! Used to compute adaptative timestep +LOGICAL :: ending_ts ! Condition to end while loop + +!---------------------------------- ! TESTS -INTEGER :: countcells -LOGICAL, SAVE :: test_flag ! flag for test/debuging outputs +INTEGER :: countcells +LOGICAL, SAVE :: test_flag ! Flag for test/debuging outputs !$OMP THREADPRIVATE(test_flag) -REAL, dimension(ngrid,nlay) :: satubf,satuaf, res_out +REAL, dimension(ngrid,nlay) :: satubf, satuaf, res_out !------------------------------------------------------------------ @@ -133,63 +133,61 @@ !============================================================= ! rad_cld is the primary radius grid used for microphysics computation. -! The grid spacing is computed assuming a constant volume ratio -! between two consecutive bins; i.e. vrat_cld. -! vrat_cld is determined from the boundary values of the size grid: -! rmin_cld and rmax_cld. +! The grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld. +! vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld. ! The rb_cld array contains the boundary values of each rad_cld bin. ! dr_cld is the width of each rad_cld bin. - ! Volume ratio between two adjacent bins - ! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3. - ! vrat_cld = exp(vrat_cld) - vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3. - vrat_cld = exp(vrat_cld) - write(*,*) "vrat_cld", vrat_cld - - rb_cld(1) = rbmin_cld - rad_cld(1) = rmin_cld - vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld - - do i=1,nbin_cld-1 - rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.) - vol_cld(i+1) = vol_cld(i) * vrat_cld - enddo - - do i=1,nbin_cld - rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cld(i) - dr_cld(i) = rb_cld(i+1) - rb_cld(i) - enddo - rb_cld(nbin_cld+1) = rbmax_cld - dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld) - - print*, ' ' - print*,'Microphysics: size bin information:' - print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)' - print*,'-----------------------------------' - do i=1,nbin_cld - write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i), dr_cld(i) - enddo - write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1) - print*,'-----------------------------------' - - ! we save that so that it is not computed at each timestep and gridpoint - rb_cld(:) = log(rb_cld(:)) - - ! Contact parameter of water ice on dust ( m=cos(theta) ) - ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - ! mteta is initialized in conf_phys - write(*,*) 'water_param contact parameter:', mteta - - ! Volume of a water molecule (m3) - vo1 = mh2o / dble(rho_ice) - ! Variance of the ice and CCN distributions - sigma_ice = sqrt(log(1.+nuice_sed)) - - write(*,*) 'Variance of ice & CCN distribs :', sigma_ice - write(*,*) 'nuice for sedimentation:', nuice_sed - write(*,*) 'Volume of a water molecule:', vo1 - - test_flag = .false. - firstcall=.false. + ! Volume ratio between two adjacent bins + ! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3. + ! vrat_cld = exp(vrat_cld) + vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3. + vrat_cld = exp(vrat_cld) + write(*,*) "vrat_cld", vrat_cld + + rb_cld(1) = rbmin_cld + rad_cld(1) = rmin_cld + vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld + + do i=1,nbin_cld-1 + rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.) + vol_cld(i+1) = vol_cld(i) * vrat_cld + enddo + + do i=1,nbin_cld + rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cld(i) + dr_cld(i) = rb_cld(i+1) - rb_cld(i) + enddo + rb_cld(nbin_cld+1) = rbmax_cld + dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld) + + write(*,*) ' ' + write(*,*)'Microphysics: size bin information:' + write(*,*)'i,rb_cld(i), rad_cld(i),dr_cld(i)' + write(*,*)'-----------------------------------' + do i=1,nbin_cld + write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i), dr_cld(i) + enddo + write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1) + write(*,*)'-----------------------------------' + + ! we save that so that it is not computed at each timestep and gridpoint + rb_cld(:) = log(rb_cld(:)) + + ! Contact parameter of water ice on dust ( m=cos(theta) ) + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! mteta is initialized in conf_phys + write(*,*) 'water_param contact parameter:', mteta + + ! Volume of a water molecule (m3) + vo1 = mh2o / dble(rho_ice) + ! Variance of the ice and CCN distributions + sigma_ice = sqrt(log(1.+nuice_sed)) + + write(*,*) 'Variance of ice & CCN distribs:', sigma_ice + write(*,*) 'nuice for sedimentation :', nuice_sed + write(*,*) 'Volume of a water molecule :', vo1 + + test_flag = .false. + firstcall = .false. ENDIF @@ -197,5 +195,4 @@ ! 1. Initialisation !============================================================= - res_out(:,:) = 0 rice(:,:) = 1.e-8 @@ -235,5 +232,4 @@ ! 2. Compute saturation !============================================================= - dev2 = 1. / ( sqrt(2.) * sigma_ice ) @@ -248,70 +244,68 @@ ! Main loop over the GCM's grid DO l=1,nlay - DO ig=1,ngrid - ! Subtimestep : here we go - IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l), zimicro(ig,l)) - spenttime = 0. - dMicetot= 0. - ending_ts=.false. - DO while (.not.ending_ts) - call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp) - zqsat(ig,l)=zqsat_tmp(1) - ! Get the partial pressure of water vapor and its saturation ratio - ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l) - satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l) - microtimestep=ptimestep/real(zimicro(ig,l)) - - ! Initialize tracers for scavenging + hdo computations (JN) - zq0(ig,l,:)=zq(ig,l,:) - - ! Check if we are integrating over ptimestep - if (spenttime+microtimestep.ge.ptimestep) then - microtimestep=ptimestep-spenttime - ! If so : last call ! - ending_ts=.true. - endif! (spenttime+microtimestep.ge.ptimestep) then + DO ig=1,ngrid + ! Subtimestep: here we go + IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l),zimicro(ig,l)) + spenttime = 0. + dMicetot= 0. + ending_ts=.false. + DO while (.not.ending_ts) + call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp) + zqsat(ig,l)=zqsat_tmp(1) + ! Get the partial pressure of water vapor and its saturation ratio + ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l) + satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l) + microtimestep=ptimestep/real(zimicro(ig,l)) + + ! Initialize tracers for scavenging + hdo computations (JN) + zq0(ig,l,:) = zq(ig,l,:) + + ! Check if we are integrating over ptimestep + if (spenttime+microtimestep >= ptimestep) then + microtimestep = ptimestep-spenttime + ! If so: last call ! + ending_ts = .true. + endif! (spenttime+microtimestep.ge.ptimestep) then !============================================================= ! 3. Nucleation !============================================================= - - IF ( satu .ge. 1. ) THEN ! if there is condensation - call updaterccn(zq(ig,l,igcm_dust_mass),zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) - - ! Expand the dust moments into a binned distribution - Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30 - No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30 - Rn = rdust(ig,l) - Rn = -log(Rn) - Rm = Rn - 3. * sigma_ice*sigma_ice - n_derf = derf( (rb_cld(1)+Rn) *dev2) - m_derf = derf( (rb_cld(1)+Rm) *dev2) - do i = 1, nbin_cld - n_aer(i) = -0.5 * No * n_derf !! this ith previously computed - m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed - n_derf = derf( (rb_cld(i+1)+Rn) *dev2) - m_derf = derf( (rb_cld(i+1)+Rm) *dev2) - n_aer(i) = n_aer(i) + 0.5 * No * n_derf - m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf - enddo - - ! Get the rates of nucleation - call nuclea(ph2o,zt(ig,l),satu,n_aer,rate) - - dN = 0. - dM = 0. - do i = 1, nbin_cld - dN = dN + n_aer(i)*(exp(-rate(i)*microtimestep)-1.) - dM = dM + m_aer(i)*(exp(-rate(i)*microtimestep)-1.) - enddo - - ! Update Dust particles - zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) - zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10) - ! Update CCNs - zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) - zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10) - - ENDIF ! of is satu >1 + IF ( satu >= 1. ) THEN ! if there is condensation + call updaterccn(zq(ig,l,igcm_dust_mass),zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) + + ! Expand the dust moments into a binned distribution + Mo = zq(ig,l,igcm_dust_mass) *tauscaling(ig) + 1.e-30 + No = zq(ig,l,igcm_dust_number)*tauscaling(ig) + 1.e-30 + Rn = rdust(ig,l) + Rn = -log(Rn) + Rm = Rn - 3. * sigma_ice*sigma_ice + n_derf = derf( (rb_cld(1)+Rn) *dev2) + m_derf = derf( (rb_cld(1)+Rm) *dev2) + do i = 1, nbin_cld + n_aer(i) = -0.5 * No * n_derf ! this ith previously computed + m_aer(i) = -0.5 * Mo * m_derf ! this ith previously computed + n_derf = derf( (rb_cld(i+1)+Rn) *dev2) + m_derf = derf( (rb_cld(i+1)+Rm) *dev2) + n_aer(i) = n_aer(i) + 0.5 * No * n_derf + m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf + enddo + + ! Get the rates of nucleation + call nuclea(ph2o,zt(ig,l),satu,n_aer,rate) + + dN = 0. + dM = 0. + do i = 1, nbin_cld + dN = dN + n_aer(i)*(exp(-rate(i)*microtimestep)-1.) + dM = dM + m_aer(i)*(exp(-rate(i)*microtimestep)-1.) + enddo + + ! Update Dust particles + zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) + zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10) + ! Update CCNs + zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) + zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10) + ENDIF ! of is satu >1 !============================================================= @@ -321,138 +315,137 @@ ! Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait ! to avoid unrealistic value for nuclei radius and so on for cases that remain negligible. - IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth - call updaterice_micro(zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),tauscaling(ig),rice(ig,l),rhocloud(ig,l)) - - No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30 - - ! saturation at equilibrium - ! rice should not be too small, otherwise seq value is not valid - seq = exp(2.*sig_h2o(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*max(rice(ig,l),1.e-7))) - - ! get resistance growth - call growthrate(zt(ig,l),pplay(ig,l),real(ph2o/satu),rice(ig,l),res,Dv) - - res_out(ig,l) = res - - ! implicit scheme of mass growth - ! cste here must be computed at each step - cste = 4*pi*rho_ice*microtimestep - - dMice = (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))/(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.) - - ! With the above scheme, dMice cannot be bigger than vapor, but can be bigger than all available ice. - dMice = max(dMice,-zq(ig,l,igcm_h2o_ice)) - dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) - - zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice - zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice - - countcells = countcells + 1 - - ! latent heat release - lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3 - subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep - - ! DIFF: trend is enforce in a range, stabilize the scheme ? - if (subpdtcloud(ig,l)*microtimestep.gt.5.0) then - subpdtcloud(ig,l)=5./microtimestep - endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then - if (subpdtcloud(ig,l)*microtimestep.lt.-5.0) then - subpdtcloud(ig,l)=-5./microtimestep - endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then - - ! Special case of the isotope of water HDO - if (hdo) then - ! condensation - if (dMice.gt.0.0) then - ! do we use fractionation? - if (hdofrac) then - ! Calculation of the HDO vapor coefficient - Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) ) * kbz * zt(ig,l) / ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo) * sqrt(1.+mhdo/mco2) ) - ! Calculation of the fractionnation coefficient at equilibrium - alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2) - ! Calculation of the 'real' fractionnation coefficient - alpha_c(ig,l) = (alpha(ig,l)*satu)/( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.) - else - alpha_c(ig,l) = 1.d0 - endif - if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then - dMice_hdo = dMice*alpha_c(ig,l)*( zq0(ig,l,igcm_hdo_vap)/zq0(ig,l,igcm_h2o_vap) ) - else - dMice_hdo=0. - endif - !! sublimation - else - if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then - dMice_hdo=dMice*( zq0(ig,l,igcm_hdo_ice)/zq0(ig,l,igcm_h2o_ice) ) - else - dMice_hdo=0. - endif - endif !if (dMice.gt.0.0) - - dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice)) - dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap)) - - zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo - zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo - - endif ! if (hdo) + IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig) >= 1.) THEN ! we trigger crystal growth + call updaterice_micro(zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),tauscaling(ig),rice(ig,l),rhocloud(ig,l)) + + No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30 + + ! saturation at equilibrium + ! rice should not be too small, otherwise seq value is not valid + seq = exp(2.*sig_h2o(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*max(rice(ig,l),1.e-7))) + + ! get resistance growth + call growthrate(zt(ig,l),pplay(ig,l),real(ph2o/satu),rice(ig,l),res,Dv) + + res_out(ig,l) = res + + ! implicit scheme of mass growth + ! cste here must be computed at each step + cste = 4*pi*rho_ice*microtimestep + + dMice = (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))/(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.) + + ! With the above scheme, dMice cannot be bigger than vapor, but can be bigger than all available ice. + dMice = max(dMice,-zq(ig,l,igcm_h2o_ice)) + dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) + + zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + dMice + zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) - dMice + + countcells = countcells + 1 + + ! latent heat release + lw = (2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3 + subpdtcloud(ig,l) = dMice*lw/cpp/microtimestep + + ! DIFF: trend is enforce in a range, stabilize the scheme? + if (subpdtcloud(ig,l)*microtimestep > 5.0) then + subpdtcloud(ig,l) = 5./microtimestep + endif + if (subpdtcloud(ig,l)*microtimestep < -5.0) then + subpdtcloud(ig,l) = -5./microtimestep + endif + + ! Special case of the isotope of water HDO + if (hdo) then + ! condensation + if (dMice > 0.) then + ! do we use fractionation? + if (hdofrac) then + ! Calculation of the HDO vapor coefficient + Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) )* kbz * zt(ig,l) / (pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo) * sqrt(1.+mhdo/mco2) ) + ! Calculation of the fractionnation coefficient at equilibrium + alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2) + ! Calculation of the 'real' fractionnation coefficient + alpha_c(ig,l) = (alpha(ig,l)*satu) / ((alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.) + else + alpha_c(ig,l) = 1.d0 + endif + if (zq0(ig,l,igcm_h2o_vap) > qperemin) then + dMice_hdo = dMice*alpha_c(ig,l)*( zq0(ig,l,igcm_hdo_vap)/zq0(ig,l,igcm_h2o_vap) ) + else + dMice_hdo = 0. + endif + ! sublimation + else + if (zq0(ig,l,igcm_h2o_ice) > qperemin) then + dMice_hdo = dMice*( zq0(ig,l,igcm_hdo_ice)/zq0(ig,l,igcm_h2o_ice) ) + else + dMice_hdo = 0. + endif + endif ! if (dMice > 0.) + + dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice)) + dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap)) + + zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice) + dMice_hdo + zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) - dMice_hdo + + endif ! if (hdo) !============================================================= ! 5. Dust cores released, tendancies, latent heat, etc ... !============================================================= - - ! If all the ice particles sublimate, all the condensation - ! nuclei are released: - if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then - ! Water - zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) + zq(ig,l,igcm_h2o_ice) - zq(ig,l,igcm_h2o_ice) = 0. - if (hdo) then - zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) + zq(ig,l,igcm_hdo_ice) - zq(ig,l,igcm_hdo_ice) = 0. - endif - ! Dust particles - zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccn_mass) - zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccn_number) - ! CCNs - zq(ig,l,igcm_ccn_mass) = 0. - zq(ig,l,igcm_ccn_number) = 0. - endif - - ELSE - ! Initialization of dMice when it's not computed - dMice=0 ! no condensation/sublimation to account for - subpdtcloud(ig,l)=0 ! no condensation/sublimation to account for - ENDIF !of if Nccn>1 - - ! No more getting tendency : we increment tracers & temp directly - - ! If not activice, the tendency from latent heat release is set to zero - IF (.not.activice) subpdtcloud(ig,l)=0. - - ! Temperature change as a feedback from latent heat release - zt(ig,l) = zt(ig,l)+subpdtcloud(ig,l)*microtimestep - - ! Prevent negative tracers ! JN - WHERE (zq(ig,l,:) < 1.e-30) zq(ig,l,:) = 1.e-30 - - IF (cloud_adapt_ts) then - ! Estimation of how much is actually condensing/subliming - dMicetot=dMicetot+abs(dMice) ! total accumulated ice formation since the beginning - IF (spenttime.ne.0) then - zdq=(dMicetot/spenttime)!*(ptimestep-spenttime) - ELSE - ! Initialization for spenttime=0 - zdq=zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep) - ENDIF - zdq=abs(zdq) - call adapt_imicro(ptimestep,zdq,zimicro(ig,l)) - ENDIF! (cloud_adapt_ts) then - ! Increment time spent in here - spenttime=spenttime+microtimestep - count_micro(ig,l)=count_micro(ig,l)+1 - ENDDO ! while (.not. ending_ts) - ENDDO ! of ig loop + ! If all the ice particles sublimate, all the condensation + ! nuclei are released: + if (zq(ig,l,igcm_h2o_ice) <= 1.e-28) then + ! Water + zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) + zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_ice) = 0. + if (hdo) then + zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) + zq(ig,l,igcm_hdo_ice) + zq(ig,l,igcm_hdo_ice) = 0. + endif + ! Dust particles + zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccn_mass) + zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccn_number) + ! CCNs + zq(ig,l,igcm_ccn_mass) = 0. + zq(ig,l,igcm_ccn_number) = 0. + endif + + ELSE + ! Initialization of dMice when it's not computed + dMice = 0 ! no condensation/sublimation to account for + subpdtcloud(ig,l) = 0 ! no condensation/sublimation to account for + ENDIF ! of if Nccn>1 + + ! No more getting tendency: we increment tracers & temp directly + + ! If not activice, the tendency from latent heat release is set to zero + IF (.not.activice) subpdtcloud(ig,l)=0. + + ! Temperature change as a feedback from latent heat release + zt(ig,l) = zt(ig,l) + subpdtcloud(ig,l)*microtimestep + + ! Prevent negative tracers ! JN + WHERE (zq(ig,l,:) < 1.e-30) zq(ig,l,:) = 1.e-30 + + IF (cloud_adapt_ts) then + ! Estimation of how much is actually condensing/subliming + dMicetot=dMicetot+abs(dMice) ! total accumulated ice formation since the beginning + IF (spenttime /= 0) then + zdq = (dMicetot/spenttime)!*(ptimestep-spenttime) + ELSE + ! Initialization for spenttime=0 + zdq = zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep) + ENDIF + zdq = abs(zdq) + call adapt_imicro(ptimestep,zdq,zimicro(ig,l)) + ENDIF ! (cloud_adapt_ts) then + ! Increment time spent in here + spenttime = spenttime + microtimestep + count_micro(ig,l) = count_micro(ig,l) + 1 + ENDDO ! while (.not. ending_ts) + ENDDO ! of ig loop ENDDO ! of nlayer loop @@ -461,13 +454,13 @@ call write_output("count_micro","count_micro after microphysics","integer",count_micro(:,:)) -!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS +!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS IF (test_flag) then - DO l=1,nlay - DO ig=1,ngrid - satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l) - satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l) + DO l=1,nlay + DO ig=1,ngrid + satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l) + satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l) + ENDDO ENDDO - ENDDO - print*, 'count is ',countcells, ' i.e. ', countcells*100/(nlay*ngrid), '% for microphys computation' + write(*,*) 'count is ',countcells, ' i.e. ', countcells*100/(nlay*ngrid), '% for microphys computation' ENDIF ! endif test_flag @@ -483,22 +476,22 @@ ! efficiency. -real,intent(in) :: ptimestep ! total duration of physics (sec) -real,intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg) -real :: alpha, beta ! Coefficients for power law -real :: defstep,coef ! Default ptimestep of 7.5 mins (iphysiq=5) -integer,intent(out) :: zimicro ! number of ptimestep division - -! Default ptimestep : defstep (7.5 mins) -defstep=88775.*5./960. -coef=ptimestep/defstep +real, intent(in) :: ptimestep ! Total duration of physics (sec) +real, intent(in) :: potcond ! Condensible vapor / sublimable ice(kg/kg) +integer, intent(out) :: zimicro ! Number of ptimestep division +real :: alpha, beta ! Coefficients for power law +real :: defstep, coef ! Default ptimestep of 7.5 mins (iphysiq=5) + +! Default ptimestep: defstep (7.5 mins) +defstep = 88775.*5./960. +coef = ptimestep/defstep ! Conservative coefficients : -! alpha=1.81846504e+11 -! beta=1.54550140e+00 -alpha=1.88282793e+05 ! Latest values for high obliquity -beta=4.57764370e-01 ! Latest values for high obliquity -!alpha=1.72198978e+10 ! Present day Mars -!beta=1.88473210e+00 -zimicro=ceiling(coef*min(max(alpha*abs(potcond)**beta,5.),7000.)) -!zimicro=2*zimicro ! Prediction times two, allow to complete year at high obliquity +!alpha = 1.81846504e+11 +!beta = 1.54550140e+00 +alpha = 1.88282793e+05 ! Latest values for high obliquity +beta = 4.57764370e-01 ! Latest values for high obliquity +!alpha = 1.72198978e+10 ! Present day Mars +!beta = 1.88473210e+00 +zimicro = ceiling(coef*min(max(alpha*abs(potcond)**beta,5.),7000.)) +!zimicro = 2*zimicro ! Prediction times two, allow to complete year at high obliquity END SUBROUTINE adapt_imicro