Changeset 4146 for trunk/LMDZ.GENERIC/libf/phygeneric/rain_generic.F90
- Timestamp:
- Mar 19, 2026, 2:35:46 PM (3 weeks ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phygeneric/rain_generic.F90
r4079 r4146 8 8 use aerosol_radius, only: aerosol_radius_h2o_liquid_ice_separate ! only used for precip_scheme_generic >=2 9 9 use tracer_h 10 use comcstfi_mod, only: g, r, cpp 10 use comcstfi_mod, only: g, cppc_ref 11 use thermo_mod 11 12 use generic_tracer_index_mod, only: generic_tracer_index 13 use callkeys_mod, only: thermo_phy,metallicity,evap_prec_generic,evap_coeff_generic, & 14 precip_scheme_generic,rainthreshold_generic,cloud_sat_generic, & 15 precip_timescale_generic,Cboucher_generic 12 16 implicit none 13 17 … … 51 55 real d_q(ngrid,nlayer) ! GCS vapor increment 52 56 real d_ql(ngrid,nlayer) ! liquid GCS / ice increment 57 real qcloud(ngrid), qcloud_tmp(ngrid) 53 58 54 59 integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice fo GCS 55 56 real, save :: RCPD ! equal to cpp57 58 real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic59 !$OMP THREADPRIVATE(metallicity)60 60 61 61 ! Subroutine options 62 62 real,parameter :: seuil_neb=0.001 ! Nebulosity threshold 63 63 64 integer,save :: precip_scheme_generic ! id number for precipitaion scheme65 66 ! for simple scheme (precip_scheme_generic=1)67 real,save :: rainthreshold_generic ! Precipitation threshold in simple scheme68 69 ! for sundquist scheme (precip_scheme_generic=2-3)70 real,save :: cloud_sat_generic ! Precipitation threshold in non simple scheme71 real,save :: precip_timescale_generic ! Precipitation timescale72 73 64 ! for Boucher scheme (precip_scheme_generic=4) 74 real,save :: Cboucher_generic ! Precipitation constant in Boucher 95 scheme75 65 real,parameter :: Kboucher=1.19E8 76 66 real,save :: c1 77 67 78 !$OMP THREADPRIVATE( precip_scheme_generic,rainthreshold_generic,cloud_sat_generic,precip_timescale_generic,Cboucher_generic,c1)68 !$OMP THREADPRIVATE(c1) 79 69 80 70 integer,parameter :: ninter=5 ! only used for precip_scheme_generic >=2 81 82 logical,save :: evap_prec_generic ! Does the rain evaporate ?83 REAL,SAVE :: evap_coeff_generic ! multiplication evaporation constant. 1. gives the model of gregory et al.84 !$OMP THREADPRIVATE(evap_prec_generic,evap_coeff_generic)85 71 86 72 ! for simple scheme : precip_scheme_generic=1 … … 102 88 real tnext(ngrid,nlayer) 103 89 104 real dmass (ngrid,nlayer) ! mass of air in each grid cell90 real dmassm2(ngrid,nlayer) ! mass of air in each grid cell 105 91 real dWtot 106 92 … … 109 95 integer, save :: i_vap_generic=0 ! GCS vapour 110 96 integer, save :: i_ice_generic=0 ! GCS ice 111 !$OMP THREADPRIVATE(i_vap_generic,i_ice_generic) 112 113 LOGICAL,save :: firstcall=.true. 114 !$OMP THREADPRIVATE(firstcall) 97 !$OMP THREADPRIVATE(i_vap_generic,i_ice_generic) 115 98 116 99 ! to call only one time the ice/vap pair of a tracer … … 142 125 RLVTT_generic=constants_RLVTT_generic(iq) 143 126 144 RCPD = cpp 145 146 147 if (firstcall) then 148 149 metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic 150 call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic 151 152 i_vap_generic=igcm_generic_vap 153 i_ice_generic=igcm_generic_ice 154 155 write(*,*) "rain: i_ice_generic=",i_ice_generic 156 write(*,*) " i_vap_generic=",i_vap_generic 157 158 write(*,*) "re-evaporate precipitations?" 159 evap_prec_generic=.true. ! default value 160 call getin_p("evap_prec_generic",evap_prec_generic) 161 write(*,*) " evap_prec_generic = ",evap_prec_generic 162 163 if (evap_prec_generic) then 164 write(*,*) "multiplicative constant in reevaporation" 165 evap_coeff_generic=1. ! default value 166 call getin_p("evap_coeff_generic",evap_coeff_generic) 167 write(*,*) " evap_coeff_generic = ",evap_coeff_generic 168 end if 169 170 write(*,*) "Precipitation scheme to use?" 171 precip_scheme_generic=1 ! default value 172 call getin_p("precip_scheme_generic",precip_scheme_generic) 173 write(*,*) " precip_scheme_generic = ",precip_scheme_generic 174 175 if (precip_scheme_generic.eq.1) then 176 write(*,*) "rainthreshold_generic in simple scheme?" 177 rainthreshold_generic=0. ! default value 178 call getin_p("rainthreshold_generic",rainthreshold_generic) 179 write(*,*) " rainthreshold_generic = ",rainthreshold_generic 180 181 !else 182 ! write(*,*) "precip_scheme_generic = ", precip_scheme_generic 183 ! write(*,*) "this precip_scheme_generic has not been implemented yet," 184 ! write(*,*) "only precip_scheme_generic = 1 has been implemented" 185 186 else if (precip_scheme_generic.eq.2.or.precip_scheme_generic.eq.3) then 187 188 write(*,*) "cloud GCS saturation level in non simple scheme?" 189 cloud_sat_generic=2.6e-4 ! default value 190 call getin_p("cloud_sat_generic",cloud_sat_generic) 191 write(*,*) " cloud_sat_generic = ",cloud_sat_generic 192 193 write(*,*) "precipitation timescale in non simple scheme?" 194 precip_timescale_generic=3600. ! default value 195 call getin_p("precip_timescale_generic",precip_timescale_generic) 196 write(*,*) " precip_timescale_generic = ",precip_timescale_generic 197 198 else if (precip_scheme_generic.eq.4) then 199 200 write(*,*) "multiplicative constant in Boucher 95 precip scheme" 201 Cboucher_generic=1. ! default value 202 call getin_p("Cboucher_generic",Cboucher_generic) 203 write(*,*) " Cboucher_generic = ",Cboucher_generic 204 205 c1=1.00*1.097/rhowater*Cboucher_generic*Kboucher 206 207 endif 208 209 PRINT*, 'in rain_generic.F, ninter=', ninter 210 211 firstcall = .false. 212 213 endif ! of if (firstcall) 127 i_vap_generic=igcm_generic_vap 128 i_ice_generic=igcm_generic_ice 214 129 215 130 ! GCM -----> subroutine variables … … 220 135 q(i,k) = pq(i,k,i_vap_generic)+pdq(i,k,i_vap_generic)*ptimestep 221 136 ql(i,k) = pq(i,k,i_ice_generic)+pdq(i,k,i_ice_generic)*ptimestep 222 223 137 !q(i,k) = pq(i,k,i_vap_generic)!+pdq(i,k,i_vap_generic) 224 138 !ql(i,k) = pq(i,k,i_ice_generic)!+pdq(i,k,i_ice_generic) … … 251 165 call Psat_generic(zt(i,k),p_tmp,metallicity,psat_tmp,zqs(i,k)) ! calculates psat_tmp & zqs(i,k) 252 166 ! metallicity --- is not used here but necessary to call function Psat_generic 253 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)167 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) 254 168 call Tsat_generic(p_tmp,Tsat(i,k)) ! calculates Tsat(i,k) 255 169 … … 260 174 do k = 1, nlayer 261 175 do i = 1, ngrid 262 dmass (i,k)=(pplev(i,k)-pplev(i,k+1))/g ! mass per m²in each layer176 dmassm2(i,k)=(pplev(i,k)-pplev(i,k+1))/g ! mass per m2 in each layer 263 177 enddo 264 178 enddo … … 275 189 276 190 if(zt(i,k).gt.Tsat(i,k))then ! if temperature of the layer box is greater than Tsat 277 ! treat the case where all liquid water should boil 278 zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*dmass(i,k)/RLVTT_generic/ptimestep,precip_rate(i)) ! on évapore 279 precip_rate(i)=MAX(precip_rate(i)-zqev,0.) ! we withdraw from precip_rate the evaporated part 280 d_q(i,k)=zqev/dmass(i,k)*ptimestep ! quantité évaporée 281 d_t(i,k) = - d_q(i,k) * RLVTT_generic/RCPD 191 ! treat the case where all condensables should boil 192 zqev=MIN((zt(i,k)-Tsat(i,k))*cppd_ref*dmassm2(i,k)/RLVTT_generic/ptimestep,precip_rate(i)) ! on évapore 193 precip_rate_tmp(i)= precip_rate(i) - zqev 194 precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0) 195 qcloud(i) = precip_rate(i)/dmassm2(i,k)*ptimestep 196 qcloud_tmp(i) = precip_rate_tmp(i)/dmassm2(i,k)*ptimestep 197 d_q(i,k)=zqev/dmassm2(i,k)*ptimestep ! quantité évaporée 198 199 SELECT CASE (TRIM(thermo_phy)) 200 CASE ('thermo_uni_ideal') 201 d_t(i,k) = d_t(i,k) - d_q(i,k) * RLVTT_generic/cppd_ref 202 END SELECT 203 precip_rate(i) = precip_rate_tmp(i) 282 204 else 205 283 206 ! zqev/zqevt are the maximum amount of vapour that we are allowed to add by evaporation of rain 284 zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmass (i,k)/(ptimestep*(1.d0+dqsat(i,k)))207 zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmassm2(i,k)/(ptimestep*(1.d0+dqsat(i,k))) 285 208 !max evaporation to reach saturation implictly accounting for temperature reduction 286 209 zqevt= MAX (0.0, evap_coeff_generic*2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5 287 *sqrt(precip_rate(i))*dmass (i,k)/pplay(i,k)*zt(i,k)*R) ! BC modif here, is it R or r/(mu/1000) ?210 *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) ? 288 211 zqev = MIN (zqev, zqevt) 289 212 zqev = MAX (zqev, 0.0) ! not necessary according to previous lines … … 292 215 precip_rate_tmp(i)= precip_rate(i) - zqev 293 216 precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0) 294 295 d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmass(i,k)*ptimestep 296 d_t(i,k) = - d_q(i,k) * RLVTT_generic/RCPD 217 qcloud(i) = precip_rate(i)/dmassm2(i,k)*ptimestep 218 qcloud_tmp(i) = precip_rate_tmp(i)/dmassm2(i,k)*ptimestep 219 220 d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmassm2(i,k)*ptimestep 221 222 SELECT CASE (TRIM(thermo_phy)) 223 224 CASE ('thermo_uni_ideal') 225 d_t(i,k) = d_t(i,k) - d_q(i,k) * RLVTT_generic/cppd_ref 226 END SELECT 227 297 228 precip_rate(i) = precip_rate_tmp(i) 229 298 230 end if 299 231 #ifdef MESOSCALE 300 d_t(i,k) = d_t(i,k)+(pphi(i,k+1)-pphi(i,k))*precip_rate(i)*ptimestep/( RCPD*dmass(i,k))232 d_t(i,k) = d_t(i,k)+(pphi(i,k+1)-pphi(i,k))*precip_rate(i)*ptimestep/(cppd_ref*dmassm2(i,k)) 301 233 ! JL22. Accounts for gravitational energy of falling precipitations (probably not to be used in the GCM 302 234 ! where the counterpart is not included in the dynamics.) … … 321 253 if ((ql(i,k)/zneb(i)).gt.lconvert)then ! precipitate! 322 254 d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0) 323 precip_rate(i) = precip_rate(i) - d_ql(i,k)*dmass (i,k)/ptimestep255 precip_rate(i) = precip_rate(i) - d_ql(i,k)*dmassm2(i,k)/ptimestep 324 256 endif 325 257 endif … … 333 265 if (rneb(i,k).GT.0.0) then 334 266 zoliq(i) = ql(i,k) 335 zrho(i) = pplay(i,k) / ( zt(i,k) * R )267 zrho(i) = pplay(i,k) / ( zt(i,k) * R(i,k) ) 336 268 zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g) 337 269 zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds) … … 415 347 if (rneb(i,k).GT.0.0) then 416 348 d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep 417 precip_rate(i) = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmass (i,k)/ptimestep349 precip_rate(i) = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmassm2(i,k)/ptimestep 418 350 endif 419 351 enddo … … 449 381 reevap_precip(i,i_vap_generic)=0. 450 382 do k=1,nlayer 451 reevap_precip(i,i_vap_generic)=reevap_precip(i,i_vap_generic)+dq_rain_generic_vap(i,k,i_vap_generic)*dmass (i,k)383 reevap_precip(i,i_vap_generic)=reevap_precip(i,i_vap_generic)+dq_rain_generic_vap(i,k,i_vap_generic)*dmassm2(i,k) 452 384 enddo 453 385 enddo
Note: See TracChangeset
for help on using the changeset viewer.
