Index: trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F90 (revision 4043) +++ trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F90 (revision 4044) @@ -79,5 +79,4 @@ REAL, dimension(ngrid,nlay) :: alpha ! HDO equilibrium fractionation coefficient (Saturation=1) REAL, dimension(ngrid,nlay) :: alpha_c ! HDO real fractionation coefficient -! REAL :: sumcheck REAL*8 :: ph2o ! Water vapor partial pressure (Pa) REAL*8 :: satu ! Water vapor saturation ratio over ice @@ -93,5 +92,5 @@ REAL :: Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient -! Parameters of the size discretization used by the microphysical scheme +! 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) @@ -111,9 +110,6 @@ REAL, dimension(ngrid,nlay) :: subpdtcloud ! Temperature variation due to latent heat REAL, dimension(ngrid,nlay,nq) :: zq0 ! local initial value of tracers -! REAL, dimension(ngrid,nlay,nq) :: zt0 ! local initial value of temperature INTEGER, dimension(ngrid,nlay) :: zimicro ! Subdivision of ptimestep INTEGER, dimension(ngrid,nlay) :: count_micro ! Number of microphys calls -REAL, dimension(ngrid,nlay) :: zpotcond_inst ! condensable water at the beginning of physics (kg/kg) -REAL, dimension(ngrid,nlay) :: zpotcond_full ! condensable water with integrated tendancies (kg/kg) REAL, dimension(ngrid,nlay) :: zpotcond ! maximal condensable water (previous two) REAL, dimension(1) :: zqsat_tmp ! maximal condensable water (previous two) @@ -154,5 +150,4 @@ rad_cld(1) = rmin_cld vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld - ! vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld do i=1,nbin_cld-1 @@ -179,5 +174,4 @@ ! we save that so that it is not computed at each timestep and gridpoint - ! rb_cld(i) = log(rb_cld(i)) rb_cld(:) = log(rb_cld(:)) @@ -207,9 +201,33 @@ rice(:,:) = 1.e-8 -! Initialize the temperature, tracers -zt(:,:)=pt(:,:) -zq(:,:,:)=pq(:,:,:) +! Initialize the temperature and tracers with tendencies + +! If scavenging, add tendency for dust all-at-once +IF (scavenging) THEN + zq(:,:,igcm_dust_mass) =zq(:,:,igcm_dust_mass)+pdq(:,:,igcm_dust_mass)*ptimestep + zq(:,:,igcm_dust_number) =zq(:,:,igcm_dust_number)+pdq(:,:,igcm_dust_number)*ptimestep +ENDIF ! scavenging + +! Add tendency for ccn all-at-once +zq(:,:,igcm_ccn_mass) = zq(:,:,igcm_ccn_mass) + pdq(:,:,igcm_ccn_mass)*ptimestep +zq(:,:,igcm_ccn_number) = zq(:,:,igcm_ccn_number) + pdq(:,:,igcm_ccn_number)*ptimestep + +! Add tendency for water all-at-once +zq(:,:,igcm_h2o_ice) = zq(:,:,igcm_h2o_ice)+pdq(:,:,igcm_h2o_ice)*ptimestep +zq(:,:,igcm_h2o_vap) = zq(:,:,igcm_h2o_vap)+pdq(:,:,igcm_h2o_vap)*ptimestep + +! Add tendency for HDO (if computed) all-at-once +IF (hdo) THEN + zq(:,:,igcm_hdo_ice) = zq(:,:,igcm_hdo_ice)+pdq(:,:,igcm_hdo_ice)*ptimestep + zq(:,:,igcm_hdo_vap) = zq(:,:,igcm_hdo_vap)+pdq(:,:,igcm_hdo_vap)*ptimestep +ENDIF + +! Add tendency for temp all-at-one +zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep + +! Local temp tendency from clouds (due to latent heath release) subpdtcloud(:,:)=0 +! Handle very small values to prevent precision error WHERE( zq(:,:,:) < 1.e-30 ) zq(:,:,:) = 1.e-30 @@ -221,22 +239,6 @@ ! Compute the condensable amount of water vapor or the sublimable water ice (if negative) -! Attention ici pdt*ptimestep peut etre severe -call watersat(ngrid*nlay,max(1.,zt+pdt*ptimestep),pplay,zqsat) ! DIFF: "temp+trend" at least 1 -zpotcond_full=(zq(:,:,igcm_h2o_vap)+pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat -zpotcond_full=max(zpotcond_full,-zq(:,:,igcm_h2o_ice)) -call watersat(ngrid*nlay,zt,pplay,zqsat) -zpotcond_inst=zq(:,:,igcm_h2o_vap) - zqsat -zpotcond_inst=max(zpotcond_inst,-zq(:,:,igcm_h2o_ice)) -! zpotcond is the most strict criterion between the two -DO l=1,nlay - DO ig=1,ngrid - if (zpotcond_full(ig,l).gt.0.) then - zpotcond(ig,l)=max(zpotcond_inst(ig,l),zpotcond_full(ig,l)) - else if (zpotcond_full(ig,l).le.0.) then - zpotcond(ig,l)=min(zpotcond_inst(ig,l),zpotcond_full(ig,l)) - endif ! (zpotcond_full.gt.0.) then - ! zpotcond(ig,l)=1.e-2 - ENDDO !ig=1,ngrid -ENDDO !l=1,nlay +call watersat(ngrid*nlay,max(1.,zt),pplay,zqsat) ! Make sure "temp+tendency" at least 1 +zpotcond=zq(:,:,igcm_h2o_vap) - zqsat countcells = 0 @@ -250,5 +252,5 @@ IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l), zimicro(ig,l)) spenttime = 0. - dMicetot= 0. ! DIFF: dMicetot=new var + dMicetot= 0. ending_ts=.false. DO while (.not.ending_ts) @@ -259,8 +261,4 @@ satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l) microtimestep=ptimestep/real(zimicro(ig,l)) - ! if (satu.lt.1.0) then - ! microtimestep=ptimestep/30. - ! write(*,*) "saturation correction" !JN - ! endif ! Initialize tracers for scavenging + hdo computations (JN) @@ -298,28 +296,4 @@ enddo - ! sumcheck = 0 - ! do i = 1, nbin_cld - ! sumcheck = sumcheck + n_aer(i) - ! enddo - ! sumcheck = abs(sumcheck/No - 1) - ! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then - ! print*, "WARNING, No sumcheck PROBLEM" - ! print*, "sumcheck, No",sumcheck, No - ! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l - ! print*, "Dust binned distribution", n_aer - ! endif - ! - ! sumcheck = 0 - ! do i = 1, nbin_cld - ! sumcheck = sumcheck + m_aer(i) - ! enddo - ! sumcheck = abs(sumcheck/Mo - 1) - ! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then - ! print*, "WARNING, Mo sumcheck PROBLEM" - ! print*, "sumcheck, Mo",sumcheck, Mo - ! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l - ! print*, "Dust binned distribution", m_aer - ! endif - ! Get the rates of nucleation call nuclea(ph2o,zt(ig,l),satu,n_aer,rate) @@ -368,8 +342,6 @@ ! 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 = max(dMice,-zq(ig,l,igcm_h2o_ice) - pdq(ig,l,igcm_h2o_ice)*microtimestep) ! DIFF: take trend into account - ! dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless... - dMice = min(dMice,zq(ig,l,igcm_h2o_vap) + pdq(ig,l,igcm_h2o_vap)*microtimestep) + 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 @@ -400,8 +372,6 @@ ! Calculation of the fractionnation coefficient at equilibrium alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2) - ! alpha = exp(13525./zt(ig,l)**2.-5.59d-2) !Lamb ! Calculation of the 'real' fractionnation coefficient alpha_c(ig,l) = (alpha(ig,l)*satu)/( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.) - ! alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics else alpha_c(ig,l) = 1.d0 @@ -459,32 +429,9 @@ ! No more getting tendency : we increment tracers & temp directly - ! Increment tracers & temp for following microtimestep - ! Dust : - ! Special treatment of dust_mass / number for scavenging ? - IF (scavenging) THEN - zq(ig,l,igcm_dust_mass) =zq(ig,l,igcm_dust_mass)+pdq(ig,l,igcm_dust_mass)*microtimestep - zq(ig,l,igcm_dust_number) =zq(ig,l,igcm_dust_number)+pdq(ig,l,igcm_dust_number)*microtimestep - ELSE - zq(ig,l,igcm_dust_mass) =zq0(ig,l,igcm_dust_mass) - zq(ig,l,igcm_dust_number) =zq0(ig,l,igcm_dust_number) - ENDIF !(scavenging) THEN - zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + pdq(ig,l,igcm_ccn_mass)*microtimestep - zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)*microtimestep - - ! Water : - zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*microtimestep - zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*microtimestep - - ! HDO (if computed) : - IF (hdo) THEN - zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+pdq(ig,l,igcm_hdo_ice)*microtimestep - zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)+pdq(ig,l,igcm_hdo_vap)*microtimestep - ENDIF ! hdo - - ! Could also set subpdtcloud to 0 if not activice to make it simpler or change name of the flag + ! If not activice, the tendency from latent heat release is set to zero IF (.not.activice) subpdtcloud(ig,l)=0. - ! Temperature : - zt(ig,l) = zt(ig,l)+(pdt(ig,l)+subpdtcloud(ig,l))*microtimestep + ! Temperature change as a feedback from latent heat release + zt(ig,l) = zt(ig,l)+subpdtcloud(ig,l)*microtimestep ! Prevent negative tracers ! JN @@ -493,17 +440,8 @@ IF (cloud_adapt_ts) then ! Estimation of how much is actually condensing/subliming - dMicetot=dMicetot+abs(dMice) ! DIFF: accumulation of absolute dMice - ! dMicetot=dMicetot+dMice - ! write(*,*) "dMicetot= ", dMicetot - ! write(*,*) "we are in (l,ig):", l,ig !JN + 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 - ! Threshold with powerlaw (sanity check) - ! zdq=min(abs(zdq), - ! & abs(zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep))) zdq=abs(zdq) call adapt_imicro(ptimestep,zdq,zimicro(ig,l)) @@ -517,6 +455,4 @@ !------ Useful outputs to check how it went -call write_output("zpotcond_inst","zpotcond_inst microphysics","(kg/kg)",zpotcond_inst(:,:)) -call write_output("zpotcond_full","zpotcond_full microphysics","(kg/kg)",zpotcond_full(:,:)) call write_output("zpotcond","zpotcond microphysics","(kg/kg)",zpotcond(:,:)) call write_output("count_micro","count_micro after microphysics","integer",count_micro(:,:)) @@ -524,56 +460,14 @@ !!!!!!!!!!!!!! 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)) 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 - print*, 'count is ',countcells, ' i.e. ', countcells*100/(nlay*ngrid), '% for microphys computation' ENDIF ! endif test_flag END SUBROUTINE improvedclouds - -!============================================================= -! The so -called "phi" function is such as phi(r) - phi(r0) = t - t0 -! It is an analytical solution to the ice radius growth equation, -! with the approximation of a constant 'reduced' cunningham correction factor -! (lambda in growthrate.F) taken at radius req instead of rice -!============================================================= - -! 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 - -!======================================================================= SUBROUTINE adapt_imicro(ptimestep,potcond,zimicro) @@ -598,8 +492,10 @@ ! alpha=1.81846504e+11 ! beta=1.54550140e+00 -alpha=1.88282793e+05 ! DIFF: new power law coefficients -beta=4.57764370e-01 -zimicro=ceiling(coef*min(max(alpha*abs(potcond)**beta,50.),7000.)) -zimicro=2*zimicro ! DIFF: prefiction times two, allow to complete year +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