Index: trunk/LMDZ.GENERIC/libf/phygeneric/rain_generic.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/rain_generic.F90 (revision 4079) +++ trunk/LMDZ.GENERIC/libf/phygeneric/rain_generic.F90 (revision 4146) @@ -8,6 +8,10 @@ use aerosol_radius, only: aerosol_radius_h2o_liquid_ice_separate ! only used for precip_scheme_generic >=2 use tracer_h - use comcstfi_mod, only: g, r, cpp + use comcstfi_mod, only: g, cppc_ref + use thermo_mod use generic_tracer_index_mod, only: generic_tracer_index + use callkeys_mod, only: thermo_phy,metallicity,evap_prec_generic,evap_coeff_generic, & + precip_scheme_generic,rainthreshold_generic,cloud_sat_generic, & + precip_timescale_generic,Cboucher_generic implicit none @@ -51,36 +55,18 @@ real d_q(ngrid,nlayer) ! GCS vapor increment real d_ql(ngrid,nlayer) ! liquid GCS / ice increment + real qcloud(ngrid), qcloud_tmp(ngrid) integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice fo GCS - - real, save :: RCPD ! equal to cpp - - real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic - !$OMP THREADPRIVATE(metallicity) ! Subroutine options real,parameter :: seuil_neb=0.001 ! Nebulosity threshold - integer,save :: precip_scheme_generic ! id number for precipitaion scheme - - ! for simple scheme (precip_scheme_generic=1) - real,save :: rainthreshold_generic ! Precipitation threshold in simple scheme - - ! for sundquist scheme (precip_scheme_generic=2-3) - real,save :: cloud_sat_generic ! Precipitation threshold in non simple scheme - real,save :: precip_timescale_generic ! Precipitation timescale - ! for Boucher scheme (precip_scheme_generic=4) - real,save :: Cboucher_generic ! Precipitation constant in Boucher 95 scheme real,parameter :: Kboucher=1.19E8 real,save :: c1 - !$OMP THREADPRIVATE(precip_scheme_generic,rainthreshold_generic,cloud_sat_generic,precip_timescale_generic,Cboucher_generic,c1) + !$OMP THREADPRIVATE(c1) integer,parameter :: ninter=5 ! only used for precip_scheme_generic >=2 - - logical,save :: evap_prec_generic ! Does the rain evaporate ? - REAL,SAVE :: evap_coeff_generic ! multiplication evaporation constant. 1. gives the model of gregory et al. - !$OMP THREADPRIVATE(evap_prec_generic,evap_coeff_generic) ! for simple scheme : precip_scheme_generic=1 @@ -102,5 +88,5 @@ real tnext(ngrid,nlayer) - real dmass(ngrid,nlayer) ! mass of air in each grid cell + real dmassm2(ngrid,nlayer) ! mass of air in each grid cell real dWtot @@ -109,8 +95,5 @@ integer, save :: i_vap_generic=0 ! GCS vapour integer, save :: i_ice_generic=0 ! GCS ice - !$OMP THREADPRIVATE(i_vap_generic,i_ice_generic) - - LOGICAL,save :: firstcall=.true. - !$OMP THREADPRIVATE(firstcall) +!$OMP THREADPRIVATE(i_vap_generic,i_ice_generic) ! to call only one time the ice/vap pair of a tracer @@ -142,74 +125,6 @@ RLVTT_generic=constants_RLVTT_generic(iq) - RCPD = cpp - - - if (firstcall) then - - metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic - call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic - - i_vap_generic=igcm_generic_vap - i_ice_generic=igcm_generic_ice - - write(*,*) "rain: i_ice_generic=",i_ice_generic - write(*,*) " i_vap_generic=",i_vap_generic - - write(*,*) "re-evaporate precipitations?" - evap_prec_generic=.true. ! default value - call getin_p("evap_prec_generic",evap_prec_generic) - write(*,*) " evap_prec_generic = ",evap_prec_generic - - if (evap_prec_generic) then - write(*,*) "multiplicative constant in reevaporation" - evap_coeff_generic=1. ! default value - call getin_p("evap_coeff_generic",evap_coeff_generic) - write(*,*) " evap_coeff_generic = ",evap_coeff_generic - end if - - write(*,*) "Precipitation scheme to use?" - precip_scheme_generic=1 ! default value - call getin_p("precip_scheme_generic",precip_scheme_generic) - write(*,*) " precip_scheme_generic = ",precip_scheme_generic - - if (precip_scheme_generic.eq.1) then - write(*,*) "rainthreshold_generic in simple scheme?" - rainthreshold_generic=0. ! default value - call getin_p("rainthreshold_generic",rainthreshold_generic) - write(*,*) " rainthreshold_generic = ",rainthreshold_generic - - !else - ! write(*,*) "precip_scheme_generic = ", precip_scheme_generic - ! write(*,*) "this precip_scheme_generic has not been implemented yet," - ! write(*,*) "only precip_scheme_generic = 1 has been implemented" - - else if (precip_scheme_generic.eq.2.or.precip_scheme_generic.eq.3) then - - write(*,*) "cloud GCS saturation level in non simple scheme?" - cloud_sat_generic=2.6e-4 ! default value - call getin_p("cloud_sat_generic",cloud_sat_generic) - write(*,*) " cloud_sat_generic = ",cloud_sat_generic - - write(*,*) "precipitation timescale in non simple scheme?" - precip_timescale_generic=3600. ! default value - call getin_p("precip_timescale_generic",precip_timescale_generic) - write(*,*) " precip_timescale_generic = ",precip_timescale_generic - - else if (precip_scheme_generic.eq.4) then - - write(*,*) "multiplicative constant in Boucher 95 precip scheme" - Cboucher_generic=1. ! default value - call getin_p("Cboucher_generic",Cboucher_generic) - write(*,*) " Cboucher_generic = ",Cboucher_generic - - c1=1.00*1.097/rhowater*Cboucher_generic*Kboucher - - endif - - PRINT*, 'in rain_generic.F, ninter=', ninter - - firstcall = .false. - - endif ! of if (firstcall) + i_vap_generic=igcm_generic_vap + i_ice_generic=igcm_generic_ice ! GCM -----> subroutine variables @@ -220,5 +135,4 @@ q(i,k) = pq(i,k,i_vap_generic)+pdq(i,k,i_vap_generic)*ptimestep ql(i,k) = pq(i,k,i_ice_generic)+pdq(i,k,i_ice_generic)*ptimestep - !q(i,k) = pq(i,k,i_vap_generic)!+pdq(i,k,i_vap_generic) !ql(i,k) = pq(i,k,i_ice_generic)!+pdq(i,k,i_ice_generic) @@ -251,5 +165,5 @@ call Psat_generic(zt(i,k),p_tmp,metallicity,psat_tmp,zqs(i,k)) ! calculates psat_tmp & zqs(i,k) ! metallicity --- is not used here but necessary to call function Psat_generic - call Lcpdqsat_generic(zt(i,k),p_tmp,psat_tmp,zqs(i,k),dqsat(i,k),dlnpsat(i,k)) ! calculates dqsat(i,k) & dlnpsat(i,k) + call Lcpdqsat_generic(zt(i,k),p_tmp,cppd_ref,psat_tmp,zqs(i,k),dqsat(i,k),dlnpsat(i,k)) ! calculates dqsat(i,k) & dlnpsat(i,k) call Tsat_generic(p_tmp,Tsat(i,k)) ! calculates Tsat(i,k) @@ -260,5 +174,5 @@ do k = 1, nlayer do i = 1, ngrid - dmass(i,k)=(pplev(i,k)-pplev(i,k+1))/g ! mass per m² in each layer + dmassm2(i,k)=(pplev(i,k)-pplev(i,k+1))/g ! mass per m2 in each layer enddo enddo @@ -275,15 +189,24 @@ if(zt(i,k).gt.Tsat(i,k))then ! if temperature of the layer box is greater than Tsat - ! treat the case where all liquid water should boil - zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*dmass(i,k)/RLVTT_generic/ptimestep,precip_rate(i)) ! on évapore - precip_rate(i)=MAX(precip_rate(i)-zqev,0.) ! we withdraw from precip_rate the evaporated part - d_q(i,k)=zqev/dmass(i,k)*ptimestep ! quantité évaporée - d_t(i,k) = - d_q(i,k) * RLVTT_generic/RCPD + ! treat the case where all condensables should boil + zqev=MIN((zt(i,k)-Tsat(i,k))*cppd_ref*dmassm2(i,k)/RLVTT_generic/ptimestep,precip_rate(i)) ! on évapore + precip_rate_tmp(i)= precip_rate(i) - zqev + precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0) + qcloud(i) = precip_rate(i)/dmassm2(i,k)*ptimestep + qcloud_tmp(i) = precip_rate_tmp(i)/dmassm2(i,k)*ptimestep + d_q(i,k)=zqev/dmassm2(i,k)*ptimestep ! quantité évaporée + + SELECT CASE (TRIM(thermo_phy)) + CASE ('thermo_uni_ideal') + d_t(i,k) = d_t(i,k) - d_q(i,k) * RLVTT_generic/cppd_ref + END SELECT + precip_rate(i) = precip_rate_tmp(i) else + ! zqev/zqevt are the maximum amount of vapour that we are allowed to add by evaporation of rain - zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmass(i,k)/(ptimestep*(1.d0+dqsat(i,k))) + zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmassm2(i,k)/(ptimestep*(1.d0+dqsat(i,k))) !max evaporation to reach saturation implictly accounting for temperature reduction zqevt= MAX (0.0, evap_coeff_generic*2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5 - *sqrt(precip_rate(i))*dmass(i,k)/pplay(i,k)*zt(i,k)*R) ! BC modif here, is it R or r/(mu/1000) ? + *sqrt(precip_rate(i))*dmassm2(i,k)/pplay(i,k)*zt(i,k)*R(i,k)) ! BC modif here, is it R or r/(mu/1000) ? zqev = MIN (zqev, zqevt) zqev = MAX (zqev, 0.0) ! not necessary according to previous lines @@ -292,11 +215,20 @@ precip_rate_tmp(i)= precip_rate(i) - zqev precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0) - - d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmass(i,k)*ptimestep - d_t(i,k) = - d_q(i,k) * RLVTT_generic/RCPD + qcloud(i) = precip_rate(i)/dmassm2(i,k)*ptimestep + qcloud_tmp(i) = precip_rate_tmp(i)/dmassm2(i,k)*ptimestep + + d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmassm2(i,k)*ptimestep + + SELECT CASE (TRIM(thermo_phy)) + + CASE ('thermo_uni_ideal') + d_t(i,k) = d_t(i,k) - d_q(i,k) * RLVTT_generic/cppd_ref + END SELECT + precip_rate(i) = precip_rate_tmp(i) + end if #ifdef MESOSCALE - d_t(i,k) = d_t(i,k)+(pphi(i,k+1)-pphi(i,k))*precip_rate(i)*ptimestep/(RCPD*dmass(i,k)) + d_t(i,k) = d_t(i,k)+(pphi(i,k+1)-pphi(i,k))*precip_rate(i)*ptimestep/(cppd_ref*dmassm2(i,k)) ! JL22. Accounts for gravitational energy of falling precipitations (probably not to be used in the GCM ! where the counterpart is not included in the dynamics.) @@ -321,5 +253,5 @@ if ((ql(i,k)/zneb(i)).gt.lconvert)then ! precipitate! d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0) - precip_rate(i) = precip_rate(i) - d_ql(i,k)*dmass(i,k)/ptimestep + precip_rate(i) = precip_rate(i) - d_ql(i,k)*dmassm2(i,k)/ptimestep endif endif @@ -333,5 +265,5 @@ if (rneb(i,k).GT.0.0) then zoliq(i) = ql(i,k) - zrho(i) = pplay(i,k) / ( zt(i,k) * R ) + zrho(i) = pplay(i,k) / ( zt(i,k) * R(i,k) ) zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g) zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds) @@ -415,5 +347,5 @@ if (rneb(i,k).GT.0.0) then d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep - precip_rate(i) = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmass(i,k)/ptimestep + precip_rate(i) = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmassm2(i,k)/ptimestep endif enddo @@ -449,5 +381,5 @@ reevap_precip(i,i_vap_generic)=0. do k=1,nlayer - reevap_precip(i,i_vap_generic)=reevap_precip(i,i_vap_generic)+dq_rain_generic_vap(i,k,i_vap_generic)*dmass(i,k) + reevap_precip(i,i_vap_generic)=reevap_precip(i,i_vap_generic)+dq_rain_generic_vap(i,k,i_vap_generic)*dmassm2(i,k) enddo enddo