Index: /trunk/LMDZ.MARS/README =================================================================== --- /trunk/LMDZ.MARS/README (revision 625) +++ /trunk/LMDZ.MARS/README (revision 626) @@ -1611,2 +1611,10 @@ cleanup around the initialization of tracer names and surface values. + +== 18/04/12 == TN +>> New scheme in improvedclouds.F based on an analytical solution of the crystal radius growth equation +>> Sedimentation of clouds done in callsedim.F, as it was done before version 520 +>> Latent heat release of sublimating ground ice in vdifc.F +>> Bug corrected in lwu.F, that was a source of NaN for very thick water ice clouds +>> Bug corrected in updatereffrad.F, ice and dust radius are now updated each timestep before radiative transfer + Index: /trunk/LMDZ.MARS/deftank/callphys.def.outliers =================================================================== --- /trunk/LMDZ.MARS/deftank/callphys.def.outliers (revision 625) +++ /trunk/LMDZ.MARS/deftank/callphys.def.outliers (revision 626) @@ -94,7 +94,7 @@ # is a cold trap) caps = .true. -# WATER: water ice albedo ? STILL TBD !! +# WATER: water ice albedo ? albedo_h2o_ice = 0.4 -# WATER: water ice thermal inertia (SI) ? STILL TBD !! +# WATER: water ice thermal inertia (SI) ? inert_h2o_ice = 1000. # WATER: water ice frost minimal thickness (ie kg.m^-2, ie 0.92 mm) for albedo @@ -107,7 +107,4 @@ # distribution of ice particles ? nuice_sed = 0.5 -# WATER: water ice cloud formation coupled with sedimentation is computed -# every "imicro" times per physical timestep -imicro = 1 ## Thermospheric options (relevant if tracer=T) : Index: /trunk/LMDZ.MARS/libf/phymars/callsedim.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/callsedim.F (revision 625) +++ /trunk/LMDZ.MARS/libf/phymars/callsedim.F (revision 626) @@ -206,15 +206,15 @@ ENDIF !of if (scavenging) -! IF (water) THEN -! write(*,*) "correction for the shape of the ice particles ?" -! beta=0.75 ! default value -! call getin("ice_shape",beta) -! write(*,*) " ice_shape = ",beta -! -! write(*,*) "water_param nueff Sedimentation:", nuice_sed -! IF (activice) THEN -! write(*,*) "water_param nueff Radiative:", nuice_ref -! ENDIF -! ENDIF + IF (water) THEN + write(*,*) "correction for the shape of the ice particles ?" + beta=0.75 ! default value + call getin("ice_shape",beta) + write(*,*) " ice_shape = ",beta + + write(*,*) "water_param nueff Sedimentation:", nuice_sed + IF (activice) THEN + write(*,*) "water_param nueff Radiative:", nuice_ref + ENDIF + ENDIF firstcall=.false. @@ -383,29 +383,29 @@ enddo ! of do ir=1,nr c ----------------------------------------------------------------- -c WATER CYCLE CASE --- NOW DONE IN WATERCLOUD !! -c ----------------------------------------------------------------- -! else if (water.and.(iq.eq.igcm_h2o_ice)) then -! if (microphys) then -! ! water ice sedimentation -! call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, -! & ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud, -! & zqi(1,1,iq),wq,beta) -! else -! ! water ice sedimentation -! call newsedim(ngrid,nlay,ngrid*nlay,1, -! & ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq), -! & zqi(1,1,iq),wq,beta) -! endif ! of if (microphys) +c WATER CYCLE CASE +c ----------------------------------------------------------------- +c else if (water.and.(iq.eq.igcm_h2o_ice)) then + else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number) + & .or. (iq .eq. igcm_h2o_ice)) then + if (microphys) then + ! water ice sedimentation + call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, + & ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud, + & zqi(1,1,iq),wq,beta) + else + ! water ice sedimentation + call newsedim(ngrid,nlay,ngrid*nlay,1, + & ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq), + & zqi(1,1,iq),wq,beta) + endif ! of if (microphys) c Tendencies c ~~~~~~~~~~ -! do ig=1,ngrid -! pdqs_sed(ig,iq)=wq(ig,1)/ptimestep -! end do + do ig=1,ngrid + pdqs_sed(ig,iq)=wq(ig,1)/ptimestep + end do c ----------------------------------------------------------------- c GENERAL CASE c ----------------------------------------------------------------- -c else - else if ((iq .ne. iccn_mass) .and. (iq .ne. iccn_number) - & .and. (iq .ne. igcm_h2o_ice)) then ! because it is done in watercloud + else call newsedim(ngrid,nlay,1,1,ptimestep, & pplev,masse,epaisseur,pt,radius(iq),rho_q(iq), @@ -445,35 +445,35 @@ c Update the ice particle size "rice" c ------------------------------------- -! IF(scavenging) THEN -! DO l = 1, nlay -! DO ig=1,ngrid -! Mo = zqi(ig,l,igcm_h2o_ice) + -! & zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30 -! No = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30 -! rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice -! & +zqi(ig,l,iccn_mass)* tauscaling(ig) -! & / Mo * rho_dust -! rhocloud(ig,l) = -! & min(max(rhocloud(ig,l),rho_ice),rho_dust) -! rice(ig,l) = -! & ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.) -! if ((Mo.lt.1.e-15) .or. (No.le.50)) rice(ig,l) = 1.e-8 -c print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No -c print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l) -! -! ENDDO -! ENDDO -! ELSE -! DO l = 1, nlay -! DO ig=1,ngrid -! ccntyp = -! & 1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.) -! ccntyp = ccntyp /ccn_factor -! rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice -! & +ccntyp*(4./3.)*pi*rdust(ig,l)**3.) -! & /(ccntyp*4./3.*pi) ), rdust(ig,l)) -! ENDDO -! ENDDO -! ENDIF + IF(scavenging) THEN + DO l = 1, nlay + DO ig=1,ngrid + Mo = zqi(ig,l,igcm_h2o_ice) + + & zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30 + No = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30 + rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice + & +zqi(ig,l,iccn_mass)* tauscaling(ig) + & / Mo * rho_dust + rhocloud(ig,l) = + & min(max(rhocloud(ig,l),rho_ice),rho_dust) + rice(ig,l) = + & ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.) + if ((Mo.lt.1.e-15) .or. (No.le.1)) rice(ig,l) = 1.e-8 +! print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No +! print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l) + + ENDDO + ENDDO + ELSE + DO l = 1, nlay + DO ig=1,ngrid + ccntyp = + & 1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.) + ccntyp = ccntyp /ccn_factor + rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice + & +ccntyp*(4./3.)*pi*rdust(ig,l)**3.) + & /(ccntyp*4./3.*pi) ), rdust(ig,l)) + ENDDO + ENDDO + ENDIF RETURN Index: /trunk/LMDZ.MARS/libf/phymars/growthrate.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/growthrate.F (revision 625) +++ /trunk/LMDZ.MARS/libf/phymars/growthrate.F (revision 626) @@ -1,4 +1,3 @@ - subroutine growthrate(timestep,temp,pmid,ph2o,psat, - & seq,rcrystal,growth) + subroutine growthrate(temp,pmid,psat,seq,rcrystal,coeff1,coeff2) IMPLICIT NONE @@ -10,4 +9,6 @@ c Authors: F. Montmessin c Adapted for the LMD/GCM by J.-B. Madeleine (October 2011) +c Use of resistances in the analytical function +c instead of growth rate - T. Navarro (2012) c c======================================================================= @@ -27,12 +28,14 @@ c ---------- - REAL timestep +c Input REAL temp ! temperature in the middle of the layer (K) REAL pmid ! pressure in the middle of the layer (K) - REAL*8 ph2o ! water vapor partial pressure (Pa) REAL*8 psat ! water vapor saturation pressure (Pa) + REAL*8 seq ! Equilibrium saturation ratio REAL rcrystal ! crystal radius before condensation (m) - REAL*8 seq ! Equilibrium saturation ratio - REAL*8 growth + +c Output + REAL coeff1,coeff2 ! coefficients for the analytical relation between time and radius + c local: @@ -43,4 +46,6 @@ REAL*8 afactor,Dv,lambda ! Intermediate computations for growth rate REAL*8 Rk,Rd + + c----------------------------------------------------------------------- @@ -73,8 +78,10 @@ & ( pi * pmid * (molco2+molh2o)*(molco2+molh2o) & * sqrt(1.+mh2o/mco2) ) - + knudsen = temp / pmid * afactor / rcrystal lambda = (1.333+0.71/knudsen) / (1.+1./knudsen) - Dv = Dv / (1. + lambda * knudsen) + +c Dv is not corrected. Instead, we use below coefficients coeff1, coeff2 +c Dv = Dv / (1. + lambda * knudsen) c - Compute Rk @@ -82,8 +89,13 @@ c - Compute Rd Rd = rgp * temp *rho_ice / (Dv*psat*mh2o) - + + + coeff1 = real(Rk + Rd*(1. + lambda * knudsen)) + coeff2 = real(Rk + Rd*(1. - lambda * knudsen)) + +c Below are growth rate used for other schemes c - Compute growth=rdr/dt, then r(t+1)= sqrt(r(t)**2.+2.*growth*dt) - growth = 1. / (seq*Rk+Rd) -c growth = (ph2o/psat-seq) / (seq*Rk+Rd) +c growth = 1. / (seq*Rk+Rd) +c growth = (ph2o/psat-seq) / (seq*Rk+Rd) c rf = sqrt( max( r**2.+2.*growth*timestep , 0. ) ) c dr = rf-r Index: /trunk/LMDZ.MARS/libf/phymars/improvedclouds.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/improvedclouds.F (revision 625) +++ /trunk/LMDZ.MARS/libf/phymars/improvedclouds.F (revision 626) @@ -1,4 +1,4 @@ subroutine improvedclouds(ngrid,nlay,ptimestep, - & pplev,pplay,zlev,pt,pdt, + & pplev,pplay,pt,pdt, & pq,pdq,pdqcloud,pdtcloud, & nq,tauscaling,rdust,rice,nuice, @@ -21,8 +21,10 @@ c values which are then used by the sedimentation and advection c schemes. +c A word about the radius growth ... +c A word about nucleation and ice growth strategies ... c Authors: J.-B. Madeleine, based on the work by Franck Montmessin c (October 2011) -c T. Navarro, debug & correction (October-December 2011) +c T. Navarro, debug,correction, new scheme (October-April 2011) c A. Spiga, optimization (February 2012) c------------------------------------------------------------------ @@ -43,5 +45,4 @@ REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) - REAL zlev(ngrid,nlay+1) ! altitude at layer boundaries REAL pt(ngrid,nlay) ! temperature at the middle of the @@ -95,9 +96,10 @@ REAL*8 Mo,No REAL*8 dN,dM,newvap - REAL*8 Rn, Rm, dev2, yeah, yeahn, yeahm + REAL*8 Rn, Rm, dev2, yeah, n_derf, m_derf REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin REAL*8 rate(nbin_cld) ! nucleation rate - REAL*8 up,dwn,Ctot,gr,seq + !REAL*8 up,dwn,Ctot,gr + REAl*8 seq REAL*8 sig ! Water-ice/air surface tension (N.m) EXTERNAL sig @@ -119,28 +121,29 @@ SAVE sigma_ice + REAL tdicho, tmax, rmin, rmax, req, rdicho + REAL coeff0, coeff1, coeff2 + REAL error_out(ngridmx,nlayermx) + REAL error2d(ngridmx) + + REAL var1,var2,var3 + + REAL rccn, epsilon + c---------------------------------- -c some outputs for 1D -- TEMPORARY +c some outputs for 1D -- TESTS REAL satu_out(ngridmx,nlayermx) ! satu ratio for output REAL dN_out(ngridmx,nlayermx) ! mass variation for output REAL dM_out(ngridmx,nlayermx) ! number variation for output - REAL dM_col(ngridmx) ! total mass condensed in column - REAL dN_col(ngridmx) ! total mass condensed in column REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!) - REAL gr_out(ngridmx,nlayermx) ! for 1d output - REAL newvap_out(ngridmx,nlayermx) ! for 1d output - REAL Mdust_col(ngridmx) ! total column dust mass - REAL Ndust_col(ngridmx) ! total column dust mass - REAL Mccn_col(ngridmx) ! total column ccn mass - REAL Nccn_col(ngridmx) ! total column ccn mass - REAL dMice_col(ngridmx) ! total column ice mass change - REAL drice_col(ngridmx) ! total column ice radius change - REAL icetot(ngridmx) ! total column ice mass + REAL gr_out(ngridmx,nlayermx) ! for 1d output REAL rice_out(ngridmx,nlayermx) ! ice radius change REAL rate_out(ngridmx,nlayermx) ! nucleation rate + REAL satubf(ngridmx,nlayermx),satuaf(ngridmx,nlayermx) + REAL ccnbf(ngridmx,nlayermx),ccnaf(ngridmx,nlayermx) INTEGER count - LOGICAL output_sca ! scavenging outputs flag for tests - - output_sca = .false. + LOGICAL test_flag ! flag for test/debuging outputs + + test_flag = .false. c---------------------------------- c---------------------------------- @@ -149,7 +152,7 @@ IF (firstcall) THEN - -c Definition of the size grid -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~ +!============================================================= +! 0. Definition of the size grid +!============================================================= c rad_cld is the primary radius grid used for microphysics computation. c The grid spacing is computed assuming a constant volume ratio @@ -215,7 +218,8 @@ END IF -c----------------------------------------------------------------------- -c 1. Initialization -c----------------------------------------------------------------------- +!============================================================= +! 1. Initialisation +!============================================================= + c Initialize the tendencies pdqcloud(1:ngrid,1:nlay,1:nq)=0 @@ -260,15 +264,16 @@ zq(ig,l,igcm_h2o_ice)= & pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep - zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004 + zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),1e-30) ! FF 12/2004 zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice) enddo enddo -c------------------------------------------------------------------ -c Cloud microphysical scheme -c------------------------------------------------------------------ - - Cste = ptimestep * 4. * pi * rho_ice +!============================================================= +! 2. Compute saturation +!============================================================= + dev2 = 1. / ( sqrt(2.) * sigma_ice ) + error_out(:,:) = 0. + call watersat(ngridmx*nlayermx,zt,pplay,zqsat) @@ -283,25 +288,28 @@ ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l) satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l) - - IF ((satu .ge. 1)! ) THEN ! if we have condensation - & .or. ( zq(ig,l,igcm_ccn_number)*tauscaling(ig) - & .ge. 2) ) THEN ! or sublimation - + + IF (( satu .ge. 1. ) .or. ! if there is condensation + & ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.)) THEN ! or sublimation + + +!============================================================= +! 3. Nucleation +!============================================================= c Expand the dust moments into a binned distribution - Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) - No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30 + 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 = -dlog(Rn) Rm = Rn - 3. * sigma_ice*sigma_ice - yeahn = derf( (rb_cld(1)+Rn) *dev2) - yeahm = derf( (rb_cld(1)+Rm) *dev2) + 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 * yeahn !! this ith previously computed - m_aer(i) = -0.5 * Mo * yeahm !! this ith previously computed - yeahn = derf( (rb_cld(i+1)+Rn) *dev2) - yeahm = derf( (rb_cld(i+1)+Rm) *dev2) - n_aer(i) = n_aer(i) + 0.5 * No * yeahn - m_aer(i) = m_aer(i) + 0.5 * Mo * yeahm + 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 !!! MORE EFFICIENT COMPUTATIONALLY THAN @@ -350,5 +358,4 @@ rate_out(ig,l) = 0 do i = 1, nbin_cld - ! schema simple n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep ) m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep ) @@ -358,85 +365,244 @@ enddo - Mo = zq0(ig,l,igcm_h2o_ice) + - & zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30 - No = zq0(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30 - !write(*,*) "l,cloud particles,cloud mass",l, No, Mo - rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice - & +zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) - & / Mo * rho_dust - rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) - if ((Mo.lt.1.e-20) .or. (No.le.1)) then - rice(ig,l) = 1.e-8 +c Update CCNs, can also be done after the radius growth +c Dust particles + zq(ig,l,igcm_dust_mass) = + & zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) + zq(ig,l,igcm_dust_number) = + & zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) +c CCNs + zq(ig,l,igcm_ccn_mass) = + & zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) + zq(ig,l,igcm_ccn_number) = + & zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) + + +!============================================================= +! 4. Ice growth: scheme for radius evolution +!============================================================= + + Mo = zq(ig,l,igcm_h2o_ice) + + & zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30 + No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30 + + rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice + & + zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + & / Mo * rho_dust + rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) + +c nuclei radius + rccn = CBRT(zq(ig,l,igcm_ccn_mass)/ + & (pi*rho_dust*zq(ig,l,igcm_ccn_number)*4./3.)) + +c ice crystal radius + rice (ig,l) = + & CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) + +c enforce physical value : crystal cannot be smaller than its nuclei ! + rice(ig,l) = max(rice(ig,l), rccn) + +c saturation at equilibrium + seq = exp( 2.*sig(zt(ig,l))*mh2o / + & (rho_ice*rgp*zt(ig,l)*rice(ig,l)) ) + + +c If there is more than on nuclei, we peform ice growth + var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN + IF (var1 .ge. -1) THEN + + + if (test_flag) then + print*, ' ' + print*, ptimestep + print*, 'satu,seq', real(satu), real(seq), ig,l + print*, 'dN,dM', real(dN),real(dM) + print*,'rccn', rccn + print*, 'Nccn, Mccn', zq(ig,l,igcm_ccn_number)*tauscaling(ig), + & zq(ig,l,igcm_ccn_mass)*tauscaling(ig) + endif + +c crystal radius to reach saturation at equilibrium (i.e. satu=seq) + req = ( zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap) + & + zq(ig,l,igcm_ccn_mass)*tauscaling(ig)* + & rho_ice/rho_dust - seq * zqsat(ig,l)) + & / ( zq(ig,l,igcm_ccn_number)*tauscaling(ig)* + & pi*rho_ice*4./3. ) + req = CBRT(req) + +c compute ice radius growth resistances (diffusive and latent heat resistancea) + call growthrate(zt(ig,l),pplay(ig,l), + & ph2o/satu,seq,req,coeff1,coeff2) + + coeff0 = -zqsat(ig,l) / (4.*pi*req*rho_ice + & * zq(ig,l,igcm_ccn_number)*tauscaling(ig)) + +c compute tmax, the time needed to reach req + call phi(rice(ig,l),req,coeff1,coeff2,tmax) + + if (test_flag) then + print*, 'coeffs', coeff0,coeff1,coeff2 + print*, 'req,tmax', req,tmax*coeff0 + print*, 'i,rmin,rdicho,rmax,tdicho' + endif + +c rmin is rice if r increases (satu >1) or req if it decreases (satu<1) +c if req is lower than rccn (ie not enough ice to reach saturation), rmin is forced to rccn + if (satu .ge. seq) then + ! crystal size is increasing + rmin = max(min(rice(ig,l),req),rccn) + rmax = max(rice(ig,l),req) else - rice(ig,l) = - & CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.) + ! crystal size is decreasing + rmax = max(min(rice(ig,l),req),rccn) + rmin = max(rice(ig,l),req) endif -c nuice(ig,l)=nuice_ref ! used for rad. transfer calculations - seq = exp( 2.*sig(zt(ig,l))*mh2o / - & (rho_ice*rgp*zt(ig,l)*rice(ig,l)) ) - - call growthrate(ptimestep,zt(ig,l),pplay(ig,l), - & ph2o,ph2o/satu,seq,rice(ig,l),gr) - - up = Cste * gr * rice(ig,l) * No * seq + - & zq(ig,l,igcm_h2o_vap) - dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1. - - Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap) - - newvap = min(up/dwn,Ctot) - - gr = gr * ( newvap/zqsat(ig,l) - seq ) - - - dMice = min( max(Cste * No * rice(ig,l) * gr, - & -zq(ig,l,igcm_h2o_ice) ), - & zq(ig,l,igcm_h2o_vap) ) - -c----------- TESTS 1D output --------- - satu_out(ig,l) = satu - Mcon_out(ig,l) = dMice - newvap_out(ig,l) = newvap - gr_out(ig,l) = gr - dN_out(ig,l) = dN - dM_out(ig,l) = dM -c------------------------------------- - -c Water ice - zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) + - & dMice - -c Water vapor - zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) - - & dMice - + !rmax = min(rmax,1.e-3) ! au max on a des rayons de 1 mm pour la dicho ... + rdicho = 0.5*(rmin+rmax) + + ! for output + var1 = rice(ig,l) + var2 = rmin + var3 = rmax + +c Given the phi function is monotonous, we perform a simple dichotomy to find the raidus at t+1 + do i = 1,10 ! dichotomy loop + +c compute tdicho, the time needed to reach rdicho + call phi(rdicho,req,coeff1,coeff2,tdicho) + !print*, tdicho,tmax + tdicho = coeff0*(tdicho - tmax) + + if (test_flag) print*, i,rmin,rdicho,rmax,tdicho + + if (tdicho .ge. ptimestep) then + rmax = rdicho + else + rmin = rdicho + endif + + rdicho = 0.5*(rmin+rmax) + + enddo ! of dichotomy loop + + if (test_flag) then + if (abs(rdicho - rccn) .ge. 1e-15) then ! to avoid infinite values + epsilon = (rmax - rmin)/(2**10) + error_out(ig,l) = + & 100*(epsilon**3 +3*epsilon**2*rdicho +3*epsilon*rdicho**2) + & / (rdicho**3-rccn**3) + endif + print*, 'error masse glace %', error_out(ig,l) + print*, 'rice,ice,vap bf', + & rice(ig,l),zq0(ig,l,igcm_h2o_ice),zq0(ig,l,igcm_h2o_vap) + endif + +c If the initial condition is subsaturated and there is not enough ice available for sublimation +c to reach equilibrium, req is neagtive. Therefore, enforce physical value. + rice(ig,l) = max(rdicho,rccn) + +!!!!! Water ice mass is computed with radius at t+1 and their current number +!!!!! Nccn is at t or t+1, depending on what has been done before +! zq(ig,l,igcm_h2o_ice) = +! & (pi*rho_ice*zq(ig,l,igcm_ccn_number)*4./3. +! & * rdicho*rdicho*rdicho - +! & zq(ig,l,igcm_ccn_mass)*rho_ice/rho_dust) +! & * tauscaling(ig) + +!!!!! Water ice mass is computed with radius at t+1 and their number at t+1 +!!!!! that is dirty, but this enforces the use of Nccn at t+1, whatever is done before +!!!!! TO BE CLEANED ONE DAY + var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN + var2 = zq0(ig,l,igcm_ccn_mass)*tauscaling(ig) + dM + zq(ig,l,igcm_h2o_ice) = + & (pi*rho_ice*var1*4./3. + & * rdicho*rdicho*rdicho - + & var2*rho_ice/rho_dust) + + + !!!! enforce realistic values (if the case of growth on Nccn(t) and condensation on Nccn(t+1)) + zq(ig,l,igcm_h2o_ice) = + & min(zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap), + & zq(ig,l,igcm_h2o_ice)) + zq(ig,l,igcm_h2o_ice) = + & max(1e-30,zq(ig,l,igcm_h2o_ice)) + + zq(ig,l,igcm_h2o_vap) = + & zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap) + & - zq(ig,l,igcm_h2o_ice) + + + if (test_flag) then + print*, 'rice,ice,vap af', + & rice(ig,l),zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_h2o_vap) + print*, 'satu bf, af', + & zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l), + & zq(ig,l,igcm_h2o_vap)/zqsat(ig,l) + endif + + +!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST + if ((zq(ig,l,igcm_h2o_ice).le. -1e-8) + & .or. (zq(ig,l,igcm_h2o_vap).le. -1e-8)) then + print*, 'NEGATIVE WATER' + print*, 'ig,l', ig,l + print*, 'satu', satu + print*, 'vap, ice bf', + & zq0(ig,l,igcm_h2o_vap), zq0(ig,l,igcm_h2o_ice) + print*, 'vap, ice af', + & zq(ig,l,igcm_h2o_vap), zq(ig,l,igcm_h2o_ice) + print*, 'ccn N,M bf', + & zq0(ig,l,igcm_ccn_number), zq0(ig,l,igcm_ccn_mass) + print*, 'ccn N,M af', + & zq(ig,l,igcm_ccn_number), zq(ig,l,igcm_ccn_mass) + print*, 'tauscaling', + & tauscaling(ig) + print*, 'req,rccn,rice bf,rice af', + & req,rccn,var1,rice(ig,l) + print*, 'rmin, rmax', var2,var3 + print*, 'error_out,tdicho,ptimestep', + & error_out(ig,l),tdicho,ptimestep + print*, ' ' + endif +!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST + + + ENDIF !of if Nccn >1 + + +!============================================================= +! 5. Dust cores released, tendancies, latent heat, etc ... +!============================================================= + c If all the ice particles sublimate, all the condensation -c nuclei are released: - if (zq(ig,l,igcm_h2o_ice).le.1e-30) then -c for coherence - dM = 0 - dN = 0 - dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig) - dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig) -c Water ice particles +c nuclei are released: + if (zq(ig,l,igcm_h2o_ice).le.1e-10) then +! for coherence +! dM = 0 +! dN = 0 +! dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig) +! dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig) +c 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. c 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) + 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) c CCNs zq(ig,l,igcm_ccn_mass) = 0. zq(ig,l,igcm_ccn_number) = 0. endif - - dN = dN/ tauscaling(ig) - dM = dM/ tauscaling(ig) -c Dust particles - zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM - zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN -c CCNs - zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM - zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN + + +! dN = dN/ tauscaling(ig) +! dM = dM/ tauscaling(ig) +!c Dust particles +! zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM +! zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN +!c CCNs +! zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM +! zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN @@ -456,4 +622,6 @@ count = count +1 + + ELSE ! for coherence (rdust, rice computations etc ..) zq(ig,l,igcm_dust_mass) = zq0(ig,l,igcm_dust_mass) @@ -463,19 +631,23 @@ zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) - + ! pour les sorties de test - satu_out(ig,l) = satu - Mcon_out(ig,l) = 0 - newvap_out(ig,l) = zq(ig,l,igcm_h2o_vap) - gr_out(ig,l) = 0 - dN_out(ig,l) = 0 - dM_out(ig,l) = 0 +! satu_out(ig,l) = satu +! gr_out(ig,l) = 0 +! dN_out(ig,l) = 0 +! dM_out(ig,l) = 0 ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice) -c-----update temperature +! ccnbf(ig,l) = zq0(ig,l,igcm_ccn_number) * tauscaling(ig) +! ccnaf(ig,l) = zq(ig,l,igcm_ccn_number) * tauscaling(ig) +! +! satubf(ig,l) = zq0(ig,l,igcm_h2o_vap) / zqsat(ig,l) +! satuaf(ig,l) = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l) + + +c-----update temperature (latent heat relase) lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3 - pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp -c pdtcloud(ig,l)=pdt(ig,l)-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp + pdtcloud(ig,l)= -pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp c----- update rice & rhocloud AFTER microphysic @@ -488,5 +660,4 @@ rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) - rice_out(ig,l)=rice(ig,l) if ((Mo.lt.1.e-20) .or. (No.le.1)) then rice(ig,l) = 1.e-8 @@ -495,5 +666,4 @@ & CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.) endif - rice_out(ig,l)=rice(ig,l)-rice_out(ig,l) nuice(ig,l)=nuice_ref ! used for rad. transfer calculations @@ -510,102 +680,63 @@ rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4) - ENDDO - ENDDO - -! print*, 'SATU' -! print*, satu_out(1,:) - -c------------------------------------------------------------------ -c------------------------------------------------------------------ -c------------------------------------------------------------------ -c------------------------------------------------------------------ -c------------------------------------------------------------------ -c TESTS - - IF (output_sca) then - + ENDDO ! of ig loop + ENDDO ! of nlayer loop + + +!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS +!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS +!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS + IF (test_flag) then + + error2d(:) = 0. + DO l=1,nlay + DO ig=1,ngrid + error2d(ig) = max(abs(error_out(ig,l)),error2d(ig)) + ENDDO + ENDDO + print*, 'count is ',count, ' i.e. ', & count*100/(nlay*ngrid), '% for microphys computation' - dM_col(:) = 0 - dN_col(:) = 0 - Mdust_col(:) = 0 - Ndust_col(:) = 0 - Mccn_col(:) = 0 - Nccn_col(:) = 0 - dMice_col(:) = 0 - drice_col(:) = 0 - icetot(:) = 0 - do l=1, nlay - do ig=1,ngrid - dM_col(ig) = dM_col(ig) + - & dM_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g - dN_col(ig) = dN_col(ig) + - & dN_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g - Mdust_col(ig) = Mdust_col(ig) + - & zq(ig,l,igcm_dust_mass)*tauscaling(ig) - & *(pplev(ig,l) - pplev(ig,l+1)) / g - Ndust_col(ig) = Ndust_col(ig) + - & zq(ig,l,igcm_dust_number)*tauscaling(ig) - & *(pplev(ig,l) - pplev(ig,l+1)) / g - Mccn_col(ig) = Mccn_col(ig) + - & zq(ig,l,igcm_ccn_mass)*tauscaling(ig) - & *(pplev(ig,l) - pplev(ig,l+1)) / g - Nccn_col(ig) = Nccn_col(ig) + - & zq(ig,l,igcm_ccn_number)*tauscaling(ig) - & *(pplev(ig,l) - pplev(ig,l+1)) / g - dMice_col(ig) = dMice_col(ig) + - & Mcon_out(ig,l) - & *(pplev(ig,l) - pplev(ig,l+1)) / g - drice_col(ig) = drice_col(ig) + - & rice_out(ig,l)*zq(ig,l,igcm_h2o_ice) - & *(pplev(ig,l) - pplev(ig,l+1)) / g - icetot(ig) = icetot(ig) + - & zq(ig,l,igcm_h2o_ice) - & *(pplev(ig,l) - pplev(ig,l+1)) / g - enddo ! of do ig=1,ngrid - enddo ! of do l=1,nlay - - drice_col(ig) = drice_col(ig)/icetot(ig) - IF (ngrid.ne.1) THEN ! 3D - call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3, - & satu_out) - call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",2, - & dM_col) - call WRITEDIAGFI(ngrid,"dN_col","dN column","N",2, - & dN_col) - call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",2, - & Ndust_col) - call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",2, - & Mdust_col) - call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",2, - & Nccn_col) - call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",2, - & Mccn_col) - call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3, - & dM_out) - call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3, - & dN_out) +! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3, +! & satu_out) +! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3, +! & dM_out) +! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3, +! & dN_out) + call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2, + & error2d) +! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3, +! & zqsat) ENDIF IF (ngrid.eq.1) THEN ! 1D - - call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1, - & newvap_out) - call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1, - & gr_out) - call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1, - & rate_out) - call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1, - & dM_out) - call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1, - & dN_out) - call WRITEDIAGFI(ngrid,"dicetot","ice col var","kg/m2",0, - & dMice_col) - call WRITEDIAGFI(ngrid,"dricetot","ice col var","m",0, - & drice_col) - call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1, - & Mcon_out) + call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1, + & error_out) + call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1, + & satubf) + call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1, + & satuaf) + call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1, + & zq0(1,:,igcm_h2o_vap)) + call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1, + & zq(1,:,igcm_h2o_vap)) + call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1, + & zq0(1,:,igcm_h2o_ice)) + call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1, + & zq(1,:,igcm_h2o_ice)) + call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1, + & ccnbf) + call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1, + & ccnaf) +c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1, +c & gr_out) +c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1, +c & rate_out) +c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1, +c & dM_out) +c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1, +c & dN_out) call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1, & zqsat) @@ -620,19 +751,55 @@ call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1, & rhocloud) - call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",0, - & dM_col) - call WRITEDIAGFI(ngrid,"dN_col","dN column","N",0, - & dN_col) - call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",0, - & Ndust_col) - call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",0, - & Mdust_col) - call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",0, - & Nccn_col) - call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",0, - & Mccn_col) ENDIF - ENDIF ! endif output_sca -c------------------------------------------------------------------ + + ENDIF ! endif test_flag +!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS +!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS +!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS + return end + + + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0 +c It is an analytical solution to the ice radius growth equation, +c with the approximation of a constant 'reduced' cunningham correction factor +c (lambda in growthrate.F) taken at radius req instead of rice +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine phi(rice,req,coeff1,coeff2,time) + + implicit none + + ! inputs + real rice ! ice radius + real req ! ice radius at equilibirum + real coeff1 ! coeff for the log + real coeff2 ! coeff for the arctan + + ! output + real time + + !local + real var + + ! 1.73205 is sqrt(3) + + var = max( + & abs(rice-req) / sqrt(rice*rice + rice*req + req*req),1e-30) + + time = + & coeff1 * + & log( var ) + & + coeff2 * 1.73205 * + & atan( (2*rice+req) / (1.73205*req) ) + + return + end + + + Index: /trunk/LMDZ.MARS/libf/phymars/lwu.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/lwu.F (revision 625) +++ /trunk/LMDZ.MARS/libf/phymars/lwu.F (revision 626) @@ -27,4 +27,6 @@ c c MODIF : FF : removing the monster bug on water ice clouds 11/2010 +c +c MODIF : TN : corrected bug if very big water ice clouds 04/2012 c----------------------------------------------------------------------- @@ -201,10 +203,14 @@ c et pourquoi pas d'abord? hourdin@lmd.ens.fr - aer_t(jl,ja,jkl)=exp(-pview*aer_a(jl,ja,jkl)) +c TN 04/12 : if very big water ice clouds, aer_t is strictly rounded +c to zero in lower levels, which is a source of NaN + !aer_t(jl,ja,jkl)=exp(-pview*aer_a(jl,ja,jkl)) + aer_t(jl,ja,jkl)=max(exp(-pview*aer_a(jl,ja,jkl)),1e-30) + enddo enddo - enddo - + enddo + c---------------------------------------------------------------------- return Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 625) +++ /trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 626) @@ -194,6 +194,4 @@ REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy - REAL watercapflag(ngridmx) ! water cap flag - c Variables used by the water ice microphysical scheme: REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m) @@ -986,5 +984,5 @@ call watercloud(ngrid,nlayer,ptimestep, & pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, - & pq,pdq,zdqcloud,zdqscloud,zdtcloud, + & pq,pdq,zdqcloud,zdtcloud, & nq,tau,tauscaling,rdust,rice,nuice, & rsedcloud,rhocloud) @@ -1008,10 +1006,14 @@ & pdq(ig,l,igcm_h2o_ice)+ & zdqcloud(ig,l,igcm_h2o_ice) - if (((pq(ig,l,igcm_h2o_ice) + + ! There are negative values issues with some tracers (17/04) + if (((pq(ig,l,igcm_h2o_ice) + & pdq(ig,l,igcm_h2o_ice)*ptimestep)) .le. 0) - & then - pdq(ig,l,igcm_h2o_ice) = - & - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30 - endif + & pdq(ig,l,igcm_h2o_ice) = + & - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30 + if (((pq(ig,l,igcm_h2o_vap) + + & pdq(ig,l,igcm_h2o_vap)*ptimestep)) .le. 0) + & pdq(ig,l,igcm_h2o_vap) = + & - pq(ig,l,igcm_h2o_vap)/ptimestep + 1e-30 + IF (scavenging) THEN pdq(ig,l,igcm_ccn_mass)= @@ -1052,10 +1054,4 @@ ENDIF ! of IF (water) THEN -! Increment water ice surface tracer tendency - DO ig=1,ngrid - dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+ - & zdqscloud(ig,igcm_h2o_ice) - ENDDO - END IF ! of IF (water) @@ -1823,13 +1819,4 @@ & 'surface h2o_ice', & 'kg.m-2',2,qsurf(1,igcm_h2o_ice)) - -c if (caps) then -c do ig=1,ngridmx -c if (watercaptag(ig)) watercapflag(ig) = 1 -c enddo -c CALL WRITEDIAGFI(ngridmx,'watercaptag', -c & 'Ice water caps', -c & '',2,watercapflag) -cs endif c CALL WRITEDIAGFI(ngridmx,'albedo', c & 'albedo', Index: /trunk/LMDZ.MARS/libf/phymars/updatereffrad.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/updatereffrad.F (revision 625) +++ /trunk/LMDZ.MARS/libf/phymars/updatereffrad.F (revision 626) @@ -69,7 +69,7 @@ REAL, PARAMETER :: ccn0 = 1.3E8 - LOGICAL firstcall - DATA firstcall/.true./ - SAVE firstcall +c LOGICAL firstcall +c DATA firstcall/.true./ +c SAVE firstcall REAL CBRT @@ -83,7 +83,10 @@ c================================================================== - IF (firstcall) THEN +c IF (firstcall) THEN c At firstcall, rdust and rice are not known; therefore c they need to be computed below. + +c Correction TN 17/04: rdust and rice must be updated at all steps, +c otherwise it is a possible source of bugs c 1.1 Dust particles @@ -120,6 +123,7 @@ ENDDO ENDIF ! of if (water.AND.activice) - firstcall = .false. - ENDIF ! of if firstcall + +c firstcall = .false. +c ENDIF ! of if firstcall c================================================================== Index: /trunk/LMDZ.MARS/libf/phymars/vdifc.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 625) +++ /trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 626) @@ -39,4 +39,5 @@ #include "comgeomfi.h" #include "tracer.h" +#include "microphys.h" c @@ -99,4 +100,6 @@ LOGICAL firstcall SAVE firstcall + + REAL lw c variable added for CO2 condensation: @@ -796,4 +799,9 @@ c Starting upward calculations for water : zq(ig,1,igcm_h2o_vap)=zq1temp(ig) +c Take into account H2o latent heat in surface energy budget + lw = (2834.3-0.28*(ptsrf(ig)-To) + & -0.004*(ptsrf(ig)-To)**2)*1.e+3 + pdtsrf(ig) = pdtsrf(ig) + & + pdqsdif(ig,igcm_h2o_ice)*lw/pcapcal(ig) ENDDO ! of DO ig=1,ngrid ELSE @@ -806,5 +814,5 @@ DO ilay=2,nlay DO ig=1,ngrid - zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) + zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) ENDDO ENDDO Index: /trunk/LMDZ.MARS/libf/phymars/watercloud.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/watercloud.F (revision 625) +++ /trunk/LMDZ.MARS/libf/phymars/watercloud.F (revision 626) @@ -1,9 +1,7 @@ - SUBROUTINE watercloud(ngrid,nlay,ptimestep, + SUBROUTINE watercloud(ngrid,nlay, ptimestep, & pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, - & pq,pdq,pdqcloud,pdqscloud,pdtcloud, + & pq,pdq,pdqcloud,pdtcloud, & nq,tau,tauscaling,rdust,rice,nuice, & rsedcloud,rhocloud) -! to use 'getin' - USE ioipsl_getincom IMPLICIT NONE @@ -15,11 +13,8 @@ c - An improved microphysical scheme (see improvedclouds.F) c -c There is a time loop specific to cloud formation and sedimentation -c due to timescales smaller than the GCM integration timestep. -c c Authors: Franck Montmessin, Francois Forget, Ehouarn Millour, c J.-B. Madeleine, Thomas Navarro c -c 2004 - 2012 +c 2004 - Oct. 2011 c======================================================================= @@ -40,28 +35,27 @@ INTEGER ngrid,nlay - INTEGER nq ! nombre de traceurs + integer nq ! nombre de traceurs REAL ptimestep ! pas de temps physique (s) REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) - REAL pdpsrf(ngrid) ! tendence surf pressure + REAL pdpsrf(ngrid) ! tendance surf pressure REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K) - REAL pdt(ngrid,nlay) ! tendence temperature des autres param. + REAL pdt(ngrid,nlay) ! tendance temperature des autres param. real pq(ngrid,nlay,nq) ! traceur (kg/kg) - real pdq(ngrid,nlay,nq) ! tendence avant condensation (kg/kg.s-1) + real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1) - REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point - REAL tauscaling(ngridmx) ! Convertion factor for dust amount - real rdust(ngridmx,nlay) ! Dust geometric mean radius (m) + REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point + REAL tauscaling(ngridmx) ! Convertion factor for dust amount + real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m) c Outputs: c ------- - real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1) - real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1) - REAL pdtcloud(ngrid,nlay) ! tendence temperature due - ! a la chaleur latente + real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1) + REAL pdtcloud(ngrid,nlay) ! tendance temperature due + ! a la chaleur latente REAL rice(ngrid,nlay) ! Ice mass mean radius (m) @@ -69,42 +63,11 @@ REAL nuice(ngrid,nlay) ! Estimated effective variance ! of the size distribution - real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius - real rhocloud(ngridmx,nlay) ! Cloud density (kg.m-3) + real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius + real rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3) c local: c ------ - ! for sedimentation - REAL zq(ngridmx,nlay,nqmx) ! local value of tracers - real masse (ngridmx,nlay) ! Layer mass (kg.m-2) - real epaisseur (ngridmx,nlay) ! Layer thickness (m) - real wq(ngridmx,nlay+1) ! displaced tracer mass (kg.m-2) - - ! for ice radius computation - REAL Mo,No - REAl ccntyp - real beta ! correction for the shape of the ice particles (cf. newsedim) - save beta - - ! for time loop - INTEGER microstep ! time subsampling step variable - INTEGER imicro ! time subsampling for coupled water microphysics & sedimentation - SAVE imicro - REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation - SAVE microtimestep - - ! tendency given by clouds (inside the micro loop) - REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud - REAL subpdtcloud(ngrid,nlay) ! cf. pdtcloud - - ! global tendency (clouds+sedim+physics) - REAL subpdq(ngrid,nlay,nq) ! cf. pdqcloud - REAL subpdt(ngrid,nlay) ! cf. pdtcloud - - REAL CBRT - EXTERNAL CBRT - - - INTEGER iq,ig,l + INTEGER ig,l LOGICAL,SAVE :: firstcall=.true. @@ -129,301 +92,22 @@ write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap write(*,*) " igcm_h2o_ice=",igcm_h2o_ice - - - write(*,*) "correction for the shape of the ice particles ?" - beta=0.75 ! default value - call getin("ice_shape",beta) - write(*,*) " ice_shape = ",beta - - write(*,*) "time subsampling for microphysic ?" - imicro = 1 - call getin("imicro",imicro) - write(*,*)"imicro = ",imicro - - microtimestep = ptimestep/float(imicro) - write(*,*)"Physical timestep is",ptimestep - write(*,*)"Microphysics timestep is",microtimestep firstcall=.false. ENDIF ! of IF (firstcall) -c-----Initialization - subpdq(1:ngrid,1:nlay,1:nq) = 0 - subpdt(1:ngrid,1:nlay) = 0 - pdqscloud(1:ngrid,1:nq) = 0 - zq(1:ngrid,1:nlay,1:nq) = 0 -c-----Computing the different layer properties for clouds sedimentation - do l=1,nlay - do ig=1, ngrid - masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g - epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l) - enddo - enddo - -c------------------------------------------------------------------ -c Time subsampling for coupled microphysics and sedimentation -c------------------------------------------------------------------ - DO microstep=1,imicro - -c------------------------------------------------------------------- -c 1. Tendencies: -c------------------ - - -c------ Temperature tendency subpdt - ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt - ! If imicro=1 subpdt is the same as pdt - DO l=1,nlay - DO ig=1,ngrid - subpdt(ig,l) = subpdt(ig,l) - & + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry - ENDDO - ENDDO -c------ Traceurs tendencies subpdq -c------ At each micro timestep we add pdq in order to have a stepped entry - IF (microphys) THEN - DO l=1,nlay - DO ig=1,ngrid - subpdq(ig,l,igcm_dust_mass) = - & subpdq(ig,l,igcm_dust_mass) - & + pdq(ig,l,igcm_dust_mass) - subpdq(ig,l,igcm_dust_number) = - & subpdq(ig,l,igcm_dust_number) - & + pdq(ig,l,igcm_dust_number) - subpdq(ig,l,igcm_ccn_mass) = - & subpdq(ig,l,igcm_ccn_mass) - & + pdq(ig,l,igcm_ccn_mass) - subpdq(ig,l,igcm_ccn_number) = - & subpdq(ig,l,igcm_ccn_number) - & + pdq(ig,l,igcm_ccn_number) - ENDDO - ENDDO - ENDIF - DO l=1,nlay - DO ig=1,ngrid - subpdq(ig,l,igcm_h2o_ice) = - & subpdq(ig,l,igcm_h2o_ice) - & + pdq(ig,l,igcm_h2o_ice) - subpdq(ig,l,igcm_h2o_vap) = - & subpdq(ig,l,igcm_h2o_vap) - & + pdq(ig,l,igcm_h2o_vap) - ENDDO - ENDDO - - -c------------------------------------------------------------------- -c 2. Main call to the different cloud schemes: -c------------------------------------------------ - IF (microphys) THEN - CALL improvedclouds(ngrid,nlay,microtimestep, - & pplev,pplay,pzlev,pt,subpdt, - & pq,subpdq,subpdqcloud,subpdtcloud, +c Main call to the different cloud schemes: + IF (microphys) THEN + CALL improvedclouds(ngrid,nlay,ptimestep, + & pplev,pplay,pt,pdt, + & pq,pdq,pdqcloud,pdtcloud, & nq,tauscaling,rdust,rice,nuice, & rsedcloud,rhocloud) - ELSE - CALL simpleclouds(ngrid,nlay,microtimestep, - & pplev,pplay,pzlev,pzlay,pt,subpdt, - & pq,subpdq,subpdqcloud,subpdtcloud, + ELSE + CALL simpleclouds(ngrid,nlay,ptimestep, + & pplev,pplay,pzlev,pzlay,pt,pdt, + & pq,pdq,pdqcloud,pdtcloud, & nq,tau,rice,nuice,rsedcloud) - ENDIF - -c-------------------------------------------------------------------- -c 3. CCN & ice sedimentation: -c-------------------------------- -! Sedimentation is done here for water ice and its CCN only -! callsedim in physics is done for dust (and others species if any) - - DO l=1,nlay - DO ig=1,ngrid - subpdt(ig,l) = - & subpdt(ig,l) + subpdtcloud(ig,l) - ENDDO - ENDDO - -c------- water ice update before sedimentation (radius is done in the cloud scheme routine) - DO l=1,nlay - DO ig=1, ngrid - zq(ig,l,igcm_h2o_ice)= max(1e-30, - & pq(ig,l,igcm_h2o_ice) + (subpdq(ig,l,igcm_h2o_ice) - & + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep) -! zq(ig,l,igcm_h2o_vap)= max(1e-30, -! & pq(ig,l,igcm_h2o_vap) + (subpdq(ig,l,igcm_h2o_vap) -! & + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep) - ENDDO - ENDDO - - -c------- CCN update before sedimentation - IF (microphys) THEN - DO l=1,nlay - DO ig=1,ngrid - zq(ig,l,igcm_ccn_number)= - & pq(ig,l,igcm_ccn_number) + (subpdq(ig,l,igcm_ccn_number) - & + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep - zq(ig,l,igcm_ccn_number)= max(1e-30, - & zq(ig,l,igcm_ccn_number))!*tauscaling(ig)) ! OU pas tauscaling ? - zq(ig,l,igcm_ccn_mass)= - & pq(ig,l,igcm_ccn_mass) + (subpdq(ig,l,igcm_ccn_mass) - & + subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep - zq(ig,l,igcm_ccn_mass)=max(1e-30, - & zq(ig,l,igcm_ccn_mass))!*tauscaling(ig)) ! OU pas tauscaling ? - ENDDO - ENDDO - ENDIF - - - - IF (microphys) THEN - -c------- CCN number sedimentation - call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, - & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, - & rhocloud,zq(1,1,igcm_ccn_number),wq,beta) - do ig=1,ngrid - ! matters if one would like to know ccn surface deposition - pdqscloud(ig,igcm_ccn_number)= - & pdqscloud(ig,igcm_ccn_number) + wq(ig,1) - enddo - -c------- CCN mass sedimentation - call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, - & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, - & rhocloud,zq(1,1,igcm_ccn_mass),wq,beta) - do ig=1,ngrid - ! matters if one would like to know ccn surface deposition - pdqscloud(ig,igcm_ccn_mass)= - & pdqscloud(ig,igcm_ccn_mass) + wq(ig,1) - enddo - -c------- Water ice sedimentation -- improved microphys - call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, - & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, - & rhocloud,zq(1,1,igcm_h2o_ice),wq,beta) - ELSE - -c------- Water ice sedimentation -- simple microphys - call newsedim(ngrid,nlay,ngrid*nlay,1, - & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, - & rho_q(igcm_h2o_ice),zq(1,1,igcm_h2o_ice),wq,beta) - - ENDIF - - -c------- Surface ice tendency update - DO ig=1,ngrid - pdqscloud(ig,igcm_h2o_ice)= - & pdqscloud(ig,igcm_h2o_ice) + wq(ig,1) - ENDDO - - -c------------------------------------------------------------------- -c 5. Updating tendencies after sedimentation: -c----------------------------------------------- - - DO l=1,nlay - DO ig=1,ngrid - - subpdq(ig,l,igcm_h2o_ice) = - & (zq(ig,l,igcm_h2o_ice) - pq(ig,l,igcm_h2o_ice)) - & /microtimestep - - - subpdq(ig,l,igcm_h2o_vap)=subpdq(ig,l,igcm_h2o_vap) - & +subpdqcloud(ig,l,igcm_h2o_vap) - - ENDDO - ENDDO - - - IF (microphys) then - DO l=1,nlay - DO ig=1,ngrid - subpdq(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number) - & -pq(ig,l,igcm_ccn_number))/microtimestep - subpdq(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass) - & -pq(ig,l,igcm_ccn_mass))/microtimestep - ENDDO - ENDDO - ENDIF - - - - ENDDO ! of DO microstep=1,imicro - -c------------------------------------------------------------------- -c 6. Compute final tendencies after time loop: -c------------------------------------------------ - -c------ Whole temperature tendency pdtcloud - DO l=1,nlay - DO ig=1,ngrid - pdtcloud(ig,l) = - & subpdt(ig,l)/imicro-pdt(ig,l) - ENDDO - ENDDO -c------ Traceurs tendencies pdqcloud - DO iq=1,nq - DO l=1,nlay - DO ig=1,ngrid - pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/imicro - & - pdq(ig,l,iq) - ENDDO - ENDDO - ENDDO -c------ Traceurs surface tendencies pdqscloud - DO iq=1,nq - DO ig=1,ngrid - pdqscloud(ig,iq) = - & pdqscloud(ig,iq)/ptimestep - ENDDO - ENDDO - - - -c------Update the ice particle size "rice" for output or photochemistry. -c------It is not used again in the water cycle. - IF(scavenging) THEN - DO l=1, nlay - DO ig=1,ngrid - Mo = pq(ig,l,igcm_h2o_ice) - & + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep - & + (pq(ig,l,igcm_ccn_mass) - & + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep) - & *tauscaling(ig) + 1.e-30 - No = (pq(ig,l,igcm_ccn_number) - & + pdqcloud(ig,l,igcm_ccn_number)*ptimestep) - & *tauscaling(ig) + 1.e-30 - rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) + - & pdqcloud(ig,l,igcm_h2o_ice)*ptimestep) - & / Mo * rho_ice - & + (pq(ig,l,igcm_ccn_mass) - & + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep) - & *tauscaling(ig)/ Mo * rho_dust - rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) - if ((Mo.lt.1.e-15) .or. (No.le.50)) then - rice(ig,l) = 1.e-8 - else - !! AS: only perform computations if conditions not met - rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.) - endif - ENDDO - ENDDO - ELSE - DO l=1,nlay - DO ig=1,ngrid - ccntyp = - & 1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.) - ccntyp = ccntyp /ccn_factor - rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice) - & + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice - & +ccntyp*(4./3.)*pi*rdust(ig,l)*rdust(ig,l)*rdust(ig,l)) - & /(ccntyp*4./3.*pi) ), rdust(ig,l)) - ENDDO - ENDDO - ENDIF ! of IF(scavenging) - -!-------------------------------------------------------------- -!-------------------------------------------------------------- + ENDIF @@ -439,6 +123,5 @@ end do -c======================================================================= - + RETURN END