Changeset 2721 for trunk/LMDZ.GENERIC/libf
- Timestamp:
- Jun 21, 2022, 11:05:45 AM (3 years ago)
- Location:
- trunk/LMDZ.GENERIC/libf/phystd
- Files:
-
- 5 edited
-
callkeys_mod.F90 (modified) (1 diff)
-
generic_cloud_common_h.F90 (modified) (1 diff)
-
inifis_mod.F90 (modified) (1 diff)
-
physiq_mod.F90 (modified) (7 diffs)
-
rain_generic.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
r2698 r2721 40 40 logical,save :: sedimentation 41 41 logical,save :: generic_condensation 42 !$OMP THREADPRIVATE(varactive,varfixed,radfixed,sedimentation,generic_condensation) 42 logical,save :: generic_rain 43 !$OMP THREADPRIVATE(varactive,varfixed,radfixed,sedimentation,generic_condensation,generic_rain) 43 44 logical,save :: water,watercond,waterrain 44 45 !$OMP THREADPRIVATE(water,watercond,waterrain) -
trunk/LMDZ.GENERIC/libf/phystd/generic_cloud_common_h.F90
r2718 r2721 206 206 207 207 if (psat .gt. p) then 208 qsat = 1. 208 qsat = 1. ! is the maximum amount of vapor that a parcel can hold without condensation, in specific concentration. 209 209 else 210 210 qsat = epsi *psat/(p-(1-epsi)*psat) -
trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
r2698 r2721 929 929 call getin_p("generic_condensation",generic_condensation) 930 930 if (is_master) write(*,*)trim(rname)//": generic_condensation = ",generic_condensation 931 932 if (is_master) write(*,*)trim(rname)//": Generic Condensation of tracers ?" 933 generic_rain=.false. !default value 934 call getin_p("generic_rain",generic_rain) 935 if (is_master) write(*,*)trim(rname)//": generic_rain = ",generic_rain 936 931 937 if (is_master) write(*,*)trim(rname)//": Compute water cycle ?" 932 938 water=.false. ! default value -
trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
r2717 r2721 61 61 startphy_file, testradtimes, tlocked, & 62 62 tracer, UseTurbDiff, water, watercond, & 63 waterrain, global1d, szangle 63 waterrain, generic_rain, global1d, szangle 64 use generic_tracer_index_mod, only: generic_tracer_index 64 65 use nonoro_gwd_ran_mod, only: nonoro_gwd_ran 65 66 use check_fields_mod, only: check_physics_fields … … 263 264 real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward) 264 265 265 integer l,ig,ierr,iq,nw,isoil,iesp 266 integer l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice 267 logical one_call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer 266 268 267 269 real zls ! Solar longitude (radians). … … 303 305 real zdtmr(ngrid,nlayer) ! Mass_redistribution routine. 304 306 real zdtrain(ngrid,nlayer) ! Rain routine. 307 real zdtrain_generic(ngrid,nlayer) ! Rain_generic routine. 305 308 real dtmoist(ngrid,nlayer) ! Moistadj routine. 306 309 real dt_ekman(ngrid,noceanmx), dt_hdiff(ngrid,noceanmx) ! Slab_ocean routine. … … 315 318 real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine. 316 319 real zdqsrain(ngrid), zdqssnow(ngrid) ! Rain routine. 320 real zdqsrain_generic(ngrid,nq), zdqssnow_generic(ngrid,nq) ! Rain_generic routine. 317 321 real dqs_hyd(ngrid,nq) ! Hydrol routine. 318 322 real reevap_precip(ngrid) ! re-evaporation flux of precipitation (integrated over the atmospheric column) 323 real reevap_precip_generic(ngrid,nq) 319 324 320 325 ! For Tracers : (kg/kg_of_air/s) … … 330 335 real dqvaplscale(ngrid,nlayer) ! Largescale routine. 331 336 real dqcldlscale(ngrid,nlayer) ! Largescale routine. 332 real dqvaplscale_generic(ngrid,nlayer,nq) !condensation_generic routine. 333 real dqcldlscale_generic(ngrid,nlayer,nq) !condensation_generic routine. 337 real dqvaplscale_generic(ngrid,nlayer,nq) ! condensation_generic routine. 338 real dqcldlscale_generic(ngrid,nlayer,nq) ! condensation_generic routine. 339 real dq_rain_generic_vap(ngrid,nlayer,nq) ! rain_generic routine 340 real dq_rain_generic_cld(ngrid,nlayer,nq) ! rain_generic routine 334 341 REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine 335 342 REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine … … 1568 1575 endif !generic_condensation 1569 1576 1577 !Generic Rain 1578 1579 if (generic_rain) then 1580 1581 zdqsrain_generic(1:ngrid,1:nq) = 0.0 1582 zdqssnow_generic(1:ngrid,1:nq) = 0.0 1583 1584 call rain_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pdt,pq,pdq, & 1585 zdtrain_generic,dq_rain_generic_vap,dq_rain_generic_cld, & 1586 zdqsrain_generic,zdqssnow_generic,reevap_precip_generic,cloudfrac) 1587 1588 pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) & 1589 + dq_rain_generic_vap(1:ngrid,1:nlayer,1:nq) 1590 ! only the parts with indexes of generic vapor tracers are filled in dq_rain_generic_vap, other parts are 0. 1591 pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) & 1592 + dq_rain_generic_cld(1:ngrid,1:nlayer,1:nq) 1593 ! only the parts with indexes of generic ice(cloud) tracers are filled in dq_rain_generic_cld, other parts are 0. 1594 1595 pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain_generic(1:ngrid,1:nlayer) 1596 1597 dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq)+zdqsrain_generic(1:ngrid,1:nq) 1598 1599 ! Test energy conservation. 1600 if(enertest)then 1601 1602 call planetwide_sumval(cpp*massarea(:,:)*zdtrain_generic(:,:)/totarea_planet,dEtot) 1603 if (is_master) print*,'In rain_generic atmospheric T energy change =',dEtot,' W m-2' 1604 1605 ! Test conservationfor each generic condensable specie (GCS) tracer 1606 1607 ! Let's loop on tracers 1608 1609 do iq=1,nq 1610 1611 one_call_ice_vap_generic = .false. 1612 1613 call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,one_call_ice_vap_generic) 1614 1615 if (one_call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer 1616 1617 call planetwide_sumval(massarea(:,:)*dq_rain_generic_cld(:,:,igcm_generic_ice)/totarea_planet*RLVTT/cpp,dItot_tmp) 1618 call planetwide_sumval(cell_area(:)*zdqssnow_generic(:,igcm_generic_ice)/totarea_planet*RLVTT/cpp,dItot) 1619 dItot = dItot + dItot_tmp 1620 call planetwide_sumval(massarea(:,:)*dq_rain_generic_vap(:,:,igcm_generic_vap)*ptimestep/totarea_planet,dVtot_tmp) 1621 call planetwide_sumval(cell_area(:)*zdqsrain_generic(:,igcm_generic_ice)/totarea_planet*RLVTT/cpp,dVtot) 1622 dVtot = dVtot + dVtot_tmp 1623 dEtot = dItot + dVtot 1624 1625 if (is_master) then 1626 print*,'In rain_generic dItot =',dItot,' W m-2' 1627 print*,'In rain_generic dVtot =',dVtot,' W m-2' 1628 print*,'In rain_generic atmospheric L energy change =',dEtot,' W m-2' 1629 endif 1630 1631 call planetwide_sumval(massarea(:,:)*dq_rain_generic_vap(:,:,igcm_generic_vap)*ptimestep/totarea_planet+ & 1632 massarea(:,:)*dq_rain_generic_cld(:,:,igcm_generic_ice)*ptimestep/totarea_planet,dWtot) 1633 call planetwide_sumval((zdqsrain_generic(:,igcm_generic_ice)+zdqssnow_generic(:,igcm_generic_ice))*cell_area(:)*ptimestep/totarea_planet,dWtots) 1634 1635 if (is_master) then 1636 print*,'In rain_generic atmospheric generic tracer change =',dWtot,' kg m-2' 1637 print*,'In rain_generic surface generic tracer change =',dWtots,' kg m-2' 1638 print*,'In rain_generic non-cons factor =',dWtot+dWtots,' kg m-2' 1639 endif 1640 1641 endif 1642 1643 end do ! do iq=1,nq loop on tracers 1644 1645 endif ! end of 'enertest' 1646 1647 endif !generic_rain 1648 1570 1649 ! Sedimentation. 1571 1650 if (sedimentation) then … … 2346 2425 call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow) 2347 2426 call writediagfi(ngrid,"reevap","reevaporation of precipitation","kg m-2 s-1",2,reevap_precip) 2427 endif 2428 2429 if(generic_rain)then 2430 call writediagfi(ngrid,"rain","generic rainfall","kg m-2 s-1",2,zdqsrain_generic) 2431 call writediagfi(ngrid,"snow","generic snowfall","kg m-2 s-1",2,zdqssnow_generic) 2432 call writediagfi(ngrid,"reevap","generic reevaporation of precipitation","kg m-2 s-1",2,reevap_precip_generic) 2348 2433 endif 2349 2434 -
trunk/LMDZ.GENERIC/libf/phystd/rain_generic.F90
r2719 r2721 1 subroutine rain_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dq rain_generic,dqsrain_generic,dqssnow_generic,reevap_precip,rneb)1 subroutine rain_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dq_rain_generic_vap,dq_rain_generic_cld,dqsrain_generic,dqssnow_generic,reevap_precip,rneb) 2 2 3 3 4 4 use ioipsl_getin_p_mod, only: getin_p 5 5 use generic_cloud_common_h 6 use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater ! T_h2O_ice_clouds,rhowater for precip_scheme_generic >=2 7 use radii_mod, only: h2o_cloudrad ! precip_scheme_generic >=2 6 use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater 7 ! T_h2O_ice_clouds,rhowater are only used for precip_scheme_generic >=2 8 use radii_mod, only: h2o_cloudrad ! only used for precip_scheme_generic >=2 8 9 use tracer_h 9 10 use comcstfi_mod, only: g, r, cpp 11 use generic_tracer_index_mod, only: generic_tracer_index 10 12 implicit none 11 13 … … 36 38 real,intent(in) :: pdq(ngrid,nlayer,nq) ! input tendency on tracers 37 39 real,intent(out) :: d_t(ngrid,nlayer) ! temperature tendency (K/s) 38 real,intent(out) :: dqrain_generic(ngrid,nlayer,nq) ! tendency of generic condensable tracers precipitation (kg/kg.s-1) 39 real,intent(out) :: dqsrain_generic(ngrid) ! rain flux at the surface (kg.m-2.s-1) 40 real,intent(out) :: dqssnow_generic(ngrid) ! snow flux at the surface (kg.m-2.s-1) 41 real,intent(out) :: reevap_precip(ngrid) ! re-evaporation flux of precipitation integrated over the atmospheric column (kg.m-2.s-1) 40 real,intent(out) :: dq_rain_generic_vap(ngrid,nlayer,nq) ! tendency of GCS tracers precipitation (kg/kg.s-1) - vapor 41 real,intent(out) :: dq_rain_generic_cld(ngrid,nlayer,nq) ! tendency of GCS tracers precipitation (kg/kg.s-1) - cloud 42 real,intent(out) :: dqsrain_generic(ngrid,nq) ! rain flux at the surface (kg.m-2.s-1) 43 real,intent(out) :: dqssnow_generic(ngrid,nq) ! snow flux at the surface (kg.m-2.s-1) 44 real,intent(out) :: reevap_precip(ngrid,nq) ! re-evaporation flux of precipitation integrated over the atmospheric column (kg.m-2.s-1) 42 45 real,intent(in) :: rneb(ngrid,nlayer) ! cloud fraction 43 46 … … 48 51 real d_ql(ngrid,nlayer) ! liquid GCS / ice increment 49 52 50 integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice of generic_tracer 51 52 real, save :: RCPD 53 real, save :: metallicity ! metallicity of planet --- is not used here but necessary to call function Psat_generic 53 integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice fo GCS 54 55 real, save :: RCPD ! equal to cpp 56 57 real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic 58 !$OMP THREADPRIVATE(metallicity) 54 59 55 60 ! Subroutine options … … 72 77 !$OMP THREADPRIVATE(precip_scheme_generic,rainthreshold_generic,cloud_sat_generic,precip_timescale_generic,Cboucher_generic,c1) 73 78 74 integer,parameter :: ninter=5 75 76 logical,save :: evap_prec_generic ! Does the rain evaporate ?79 integer,parameter :: ninter=5 ! only used for precip_scheme_generic >=2 80 81 logical,save :: evap_prec_generic ! Does the rain evaporate ? 77 82 REAL,SAVE :: evap_coeff_generic ! multiplication evaporation constant. 1. gives the model of gregory et al. 78 83 !$OMP THREADPRIVATE(evap_prec_generic evap_coeff_generic) 79 84 80 ! for simple scheme 85 ! for simple scheme : precip_scheme_generic=1 81 86 real,parameter :: t_crit=218.0 ! need to be modified to be generic for all species 82 87 real lconvert … … 109 114 !$OMP THREADPRIVATE(firstcall) 110 115 116 ! to call only one time the ice/vap pair of a tracer 117 logical one_call_ice_vap_generic 118 111 119 ! Online functions 112 120 real fallv, fall2v, zzz ! falling speed of ice crystals … … 114 122 fall2v (zzz) =10.6 * ((zzz)**0.31) !for use with radii 115 123 116 ! Let's loop on tracers 124 ! Let's loop on tracers 117 125 118 126 do iq=1,nq 119 if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0) & 120 .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"co2") .eq. 0)) then 127 128 one_call_ice_vap_generic = .false. 129 130 call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,one_call_ice_vap_generic) 131 132 if(one_call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer 133 121 134 ! Let's get the index of our tracers (we look for igcm _generic_vap and igcm_generic_ice) 122 135 igcm_generic_vap=iq … … 135 148 end if 136 149 if (igcm_generic_ice .eq. -1) then 137 138 150 write(*,*) "ERROR : You set a vap traceur but you forgot to set the corresponding ice traceur, & 151 or the pair vap/ice is not written one after another in traceur.def" 139 152 endif 140 end if ! if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0))141 142 m=constants_mass(iq)143 delta_vapH=constants_delta_vapH(iq)144 Tref=constants_Tref(iq)145 Pref=constants_Pref(iq)146 epsi=constants_epsi(iq)147 RLVTT=constants_RLVTT(iq)148 149 RCPD = cpp150 153 151 end do ! do iq=1,nq loop on tracers 152 153 if (firstcall) then 154 155 metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic 156 call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic 157 158 i_vap_generic=igcm_generic_vap 159 i_ice_generic=igcm_generic_ice 160 161 write(*,*) "rain: i_ice_generic=",i_ice_generic 162 write(*,*) " i_vap_generic=",i_vap_generic 163 164 PRINT*, 'in rain_generic.F, ninter=', ninter 165 PRINT*, 'in rain_generic.F, evap_prec_generic=', evap_prec_generic 166 167 write(*,*) "Precipitation scheme to use?" 168 precip_scheme_generic=1 ! default value 169 call getin_p("precip_scheme_generic",precip_scheme_generic) 170 write(*,*) " precip_scheme_generic = ",precip_scheme_generic 171 172 if (precip_scheme_generic.eq.1) then 173 write(*,*) "rainthreshold_generic in simple scheme?" 174 rainthreshold_generic=0. ! default value 175 call getin_p("rainthreshold_generic",rainthreshold_generic) 176 write(*,*) " rainthreshold_generic = ",rainthreshold_generic 177 178 !else 179 ! write(*,*) "precip_scheme_generic = ", precip_scheme_generic 180 ! write(*,*) "this precip_scheme_generic has not been implemented yet," 181 ! write(*,*) "only precip_scheme_generic = 1 has been implemented" 154 m=constants_mass(iq) 155 delta_vapH=constants_delta_vapH(iq) 156 Tref=constants_Tref(iq) 157 Pref=constants_Pref(iq) 158 epsi=constants_epsi(iq) 159 RLVTT=constants_RLVTT(iq) 160 161 RCPD = cpp 182 162 183 else if (precip_scheme_generic.eq.2.or.precip_scheme_generic.eq.3) then 184 185 write(*,*) "cloud GCS saturation level in non simple scheme?" 186 cloud_sat_generic=2.6e-4 ! default value 187 call getin_p("cloud_sat_generic",cloud_sat_generic) 188 write(*,*) " cloud_sat_generic = ",cloud_sat_generic 189 190 write(*,*) "precipitation timescale in non simple scheme?" 191 precip_timescale_generic=3600. ! default value 192 call getin_p("precip_timescale_generic",precip_timescale_generic) 193 write(*,*) " precip_timescale_generic = ",precip_timescale_generic 194 195 else if (precip_scheme_generic.eq.4) then 196 197 write(*,*) "multiplicative constant in Boucher 95 precip scheme" 198 Cboucher_generic=1. ! default value 199 call getin_p("Cboucher_generic",Cboucher_generic) 200 write(*,*) " Cboucher_generic = ",Cboucher_generic 201 202 c1=1.00*1.097/rhowater*Cboucher_generic*Kboucher 203 204 endif 205 206 write(*,*) "re-evaporate precipitations?" 207 evap_prec_generic=.true. ! default value 208 call getin_p("evap_prec_generic",evap_prec_generic) 209 write(*,*) " evap_prec_generic = ",evap_prec_generic 210 211 if (evap_prec_generic) then 212 write(*,*) "multiplicative constant in reevaporation" 213 evap_coeff_generic=1. ! default value 214 call getin_p("evap_coeff_generic",evap_coeff_generic) 215 write(*,*) " evap_coeff_generic = ",evap_coeff_generic 216 end if 217 218 firstcall = .false. 219 220 endif ! of if (firstcall) 221 222 ! GCM -----> subroutine variables 223 do k = 1, nlayer 224 do i = 1, ngrid 225 226 zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here 227 q(i,k) = pq(i,k,i_vap_generic)+pdq(i,k,i_vap_generic)*ptimestep 228 ql(i,k) = pq(i,k,i_ice_generic)+pdq(i,k,i_ice_generic)*ptimestep 229 230 !q(i,k) = pq(i,k,i_vap_generic)!+pdq(i,k,i_vap_generic) 231 !ql(i,k) = pq(i,k,i_ice_generic)!+pdq(i,k,i_ice_generic) 232 233 if(q(i,k).lt.0.)then ! if this is not done, we don't conserve GCS 234 q(i,k)=0. 235 endif 236 if(ql(i,k).lt.0.)then 237 ql(i,k)=0. 238 endif 239 enddo 240 enddo 241 242 ! Initialise the outputs 243 d_t(1:ngrid,1:nlayer) = 0.0 244 d_q(1:ngrid,1:nlayer) = 0.0 245 d_ql(1:ngrid,1:nlayer) = 0.0 246 precip_rate(1:ngrid) = 0.0 247 precip_rate_tmp(1:ngrid) = 0.0 248 249 ! calculate saturation mixing ratio 250 do k = 1, nlayer 251 do i = 1, ngrid 252 p_tmp = pplay(i,k) 253 call Psat_generic(zt(i,k),p_tmp,metallicity,psat_tmp,zqs(i,k)) ! metallicit --- is not used here but necessary to call function Psat_generic 254 call Lcpdqsat_generic(zt(i,k),p_tmp,psat_tmp,zqs(i,k),dqsat(i,k),dlnpsat(i,k)) 255 call Tsat_generic(p_tmp,Tsat(i,k)) 256 257 enddo 258 enddo 259 260 ! get column / layer conversion factor 261 do k = 1, nlayer 262 do i = 1, ngrid 263 dmass(i,k)=(pplev(i,k)-pplev(i,k+1))/g 264 enddo 265 enddo 266 267 ! Vertical loop (from top to bottom) 268 ! We carry the rain with us and calculate that added by warm/cold precipitation 269 ! processes and that subtracted by evaporation at each level. 270 do k = nlayer, 1, -1 271 272 if (evap_prec_generic) then ! note no rneb dependence! 273 do i = 1, ngrid 274 if (precip_rate(i) .GT.0.) then 275 276 if(zt(i,k).gt.Tsat(i,k))then 277 ! treat the case where all liquid water should boil 278 zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*dmass(i,k)/RLVTT/ptimestep,precip_rate(i)) ! on évapore 279 precip_rate(i)=MAX(precip_rate(i)-zqev,0.) 280 d_q(i,k)=zqev/dmass(i,k)*ptimestep ! quantité évaporée 281 d_t(i,k) = - d_q(i,k) * RLVTT/RCPD 282 else 283 zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmass(i,k)/(ptimestep*(1.d0+dqsat(i,k))) 284 !max evaporation to reach saturation implictly accounting for temperature reduction 285 zqevt= MAX (0.0, evap_coeff_generic*2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5 286 *sqrt(precip_rate(i))*dmass(i,k)/pplay(i,k)*zt(i,k)*R) ! BC modif here 287 zqev = MIN (zqev, zqevt) 288 zqev = MAX (zqev, 0.0) 289 precip_rate_tmp(i)= precip_rate(i) - zqev 290 precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0) 291 292 d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmass(i,k)*ptimestep 293 d_t(i,k) = - d_q(i,k) * RLVTT/RCPD 294 precip_rate(i) = precip_rate_tmp(i) 295 end if 296 endif ! of if (precip_rate(i) .GT.0.) 163 164 if (firstcall) then 165 166 metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic 167 call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic 168 169 i_vap_generic=igcm_generic_vap 170 i_ice_generic=igcm_generic_ice 171 172 write(*,*) "rain: i_ice_generic=",i_ice_generic 173 write(*,*) " i_vap_generic=",i_vap_generic 174 175 write(*,*) "re-evaporate precipitations?" 176 evap_prec_generic=.true. ! default value 177 call getin_p("evap_prec_generic",evap_prec_generic) 178 write(*,*) " evap_prec_generic = ",evap_prec_generic 179 180 if (evap_prec_generic) then 181 write(*,*) "multiplicative constant in reevaporation" 182 evap_coeff_generic=1. ! default value 183 call getin_p("evap_coeff_generic",evap_coeff_generic) 184 write(*,*) " evap_coeff_generic = ",evap_coeff_generic 185 end if 186 187 write(*,*) "Precipitation scheme to use?" 188 precip_scheme_generic=1 ! default value 189 call getin_p("precip_scheme_generic",precip_scheme_generic) 190 write(*,*) " precip_scheme_generic = ",precip_scheme_generic 191 192 if (precip_scheme_generic.eq.1) then 193 write(*,*) "rainthreshold_generic in simple scheme?" 194 rainthreshold_generic=0. ! default value 195 call getin_p("rainthreshold_generic",rainthreshold_generic) 196 write(*,*) " rainthreshold_generic = ",rainthreshold_generic 197 198 !else 199 ! write(*,*) "precip_scheme_generic = ", precip_scheme_generic 200 ! write(*,*) "this precip_scheme_generic has not been implemented yet," 201 ! write(*,*) "only precip_scheme_generic = 1 has been implemented" 202 203 else if (precip_scheme_generic.eq.2.or.precip_scheme_generic.eq.3) then 204 205 write(*,*) "cloud GCS saturation level in non simple scheme?" 206 cloud_sat_generic=2.6e-4 ! default value 207 call getin_p("cloud_sat_generic",cloud_sat_generic) 208 write(*,*) " cloud_sat_generic = ",cloud_sat_generic 209 210 write(*,*) "precipitation timescale in non simple scheme?" 211 precip_timescale_generic=3600. ! default value 212 call getin_p("precip_timescale_generic",precip_timescale_generic) 213 write(*,*) " precip_timescale_generic = ",precip_timescale_generic 214 215 else if (precip_scheme_generic.eq.4) then 216 217 write(*,*) "multiplicative constant in Boucher 95 precip scheme" 218 Cboucher_generic=1. ! default value 219 call getin_p("Cboucher_generic",Cboucher_generic) 220 write(*,*) " Cboucher_generic = ",Cboucher_generic 221 222 c1=1.00*1.097/rhowater*Cboucher_generic*Kboucher 223 224 endif 225 226 PRINT*, 'in rain_generic.F, ninter=', ninter 227 228 firstcall = .false. 229 230 endif ! of if (firstcall) 231 232 ! GCM -----> subroutine variables 233 do k = 1, nlayer 234 do i = 1, ngrid 235 236 zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here 237 q(i,k) = pq(i,k,i_vap_generic)+pdq(i,k,i_vap_generic)*ptimestep 238 ql(i,k) = pq(i,k,i_ice_generic)+pdq(i,k,i_ice_generic)*ptimestep 239 240 !q(i,k) = pq(i,k,i_vap_generic)!+pdq(i,k,i_vap_generic) 241 !ql(i,k) = pq(i,k,i_ice_generic)!+pdq(i,k,i_ice_generic) 242 243 if(q(i,k).lt.0.)then ! if this is not done, we don't conserve GCS 244 q(i,k)=0. ! vap 245 endif 246 if(ql(i,k).lt.0.)then 247 ql(i,k)=0. ! ice 248 endif 249 enddo 297 250 enddo 298 endif ! of if (evap_prec_generic) 299 300 zoliq(1:ngrid) = 0.0 301 302 303 if(precip_scheme_generic.eq.1)then 304 305 do i = 1, ngrid 306 t_tmp = zt(i,k) 307 if (t_tmp .ge. T_h2O_ice_liq) then 308 lconvert=rainthreshold_generic 309 elseif (t_tmp .gt. t_crit) then 310 lconvert=rainthreshold_generic*(1.- t_crit/t_tmp) 311 lconvert=MAX(0.0,lconvert) 312 else 313 lconvert=0. 314 endif 315 316 317 if (ql(i,k).gt.1.e-9) then 318 zneb(i) = MAX(rneb(i,k), seuil_neb) 319 if ((ql(i,k)/zneb(i)).gt.lconvert)then ! precipitate! 320 d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0) 321 precip_rate(i) = precip_rate(i) - d_ql(i,k)*dmass(i,k)/ptimestep 322 endif 323 endif 251 252 ! Initialise the outputs 253 d_t(1:ngrid,1:nlayer) = 0.0 254 d_q(1:ngrid,1:nlayer) = 0.0 255 d_ql(1:ngrid,1:nlayer) = 0.0 256 precip_rate(1:ngrid) = 0.0 257 precip_rate_tmp(1:ngrid) = 0.0 258 259 ! calculate saturation mixing ratio 260 do k = 1, nlayer 261 do i = 1, ngrid 262 p_tmp = pplay(i,k) 263 call Psat_generic(zt(i,k),p_tmp,metallicity,psat_tmp,zqs(i,k)) ! calculates psat_tmp & zqs(i,k) 264 ! metallicity --- is not used here but necessary to call function Psat_generic 265 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) 266 call Tsat_generic(p_tmp,Tsat(i,k)) ! calculates Tsat(i,k) 267 268 enddo 324 269 enddo 325 270 326 elseif (precip_scheme_generic.ge.2) then 327 328 do i = 1, ngrid 329 if (rneb(i,k).GT.0.0) then 330 zoliq(i) = ql(i,k) 331 zrho(i) = pplay(i,k) / ( zt(i,k) * R ) 332 zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g) 333 zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds) 334 zfrac(i) = MAX(zfrac(i), 0.0) 335 zfrac(i) = MIN(zfrac(i), 1.0) 336 zneb(i) = MAX(rneb(i,k), seuil_neb) 337 endif 271 ! get column / layer conversion factor 272 do k = 1, nlayer 273 do i = 1, ngrid 274 dmass(i,k)=(pplev(i,k)-pplev(i,k+1))/g ! mass per m² in each layer 275 enddo 338 276 enddo 339 277 340 !recalculate liquid GCS particle radii 341 call h2o_cloudrad(ngrid,nlayer,ql,reffh2oliq,reffh2oice) 342 343 SELECT CASE(precip_scheme_generic) 344 345 !precip scheme from Sundquist 78 346 CASE(2) 347 348 do n = 1, ninter 278 ! Vertical loop (from top to bottom) 279 ! We carry the rain with us and calculate that added by warm/cold precipitation 280 ! processes and that subtracted by evaporation at each level. 281 ! We go from a layer to the next layer below and make the rain fall 282 do k = nlayer, 1, -1 283 284 if (evap_prec_generic) then ! note no rneb dependence! 285 do i = 1, ngrid 286 if (precip_rate(i) .GT.0.) then ! if rain from upper layers has fallen in the current layer box 287 288 if(zt(i,k).gt.Tsat(i,k))then ! if temperature of the layer box is greater than Tsat 289 ! treat the case where all liquid water should boil 290 zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*dmass(i,k)/RLVTT/ptimestep,precip_rate(i)) ! on évapore 291 precip_rate(i)=MAX(precip_rate(i)-zqev,0.) ! we withdraw from precip_rate the evaporated part 292 d_q(i,k)=zqev/dmass(i,k)*ptimestep ! quantité évaporée 293 d_t(i,k) = - d_q(i,k) * RLVTT/RCPD 294 else 295 zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmass(i,k)/(ptimestep*(1.d0+dqsat(i,k))) 296 !max evaporation to reach saturation implictly accounting for temperature reduction 297 zqevt= MAX (0.0, evap_coeff_generic*2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5 298 *sqrt(precip_rate(i))*dmass(i,k)/pplay(i,k)*zt(i,k)*R) ! BC modif here, is it R or r/(mu/1000) ? 299 zqev = MIN (zqev, zqevt) 300 zqev = MAX (zqev, 0.0) ! a priori inutile d'après les précédentes lignes 301 precip_rate_tmp(i)= precip_rate(i) - zqev 302 precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0) 303 304 d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmass(i,k)*ptimestep 305 d_t(i,k) = - d_q(i,k) * RLVTT/RCPD 306 precip_rate(i) = precip_rate_tmp(i) 307 end if 308 endif ! of if (precip_rate(i) .GT.0.) 309 enddo 310 endif ! of if (evap_prec_generic) 311 312 zoliq(1:ngrid) = 0.0 313 314 315 if(precip_scheme_generic.eq.1)then 316 317 do i = 1, ngrid 318 t_tmp = zt(i,k) 319 if (t_tmp .ge. T_h2O_ice_liq) then 320 lconvert=rainthreshold_generic 321 elseif (t_tmp .gt. t_crit) then 322 lconvert=rainthreshold_generic*(1.- t_crit/t_tmp) 323 lconvert=MAX(0.0,lconvert) 324 else 325 lconvert=0. 326 endif 327 328 ! lconvert is between 0. and rainthreshold_generic 329 330 331 if (ql(i,k).gt.1.e-9) then 332 zneb(i) = MAX(rneb(i,k), seuil_neb) ! in mesoscale rneb = 0 or 1 333 if ((ql(i,k)/zneb(i)).gt.lconvert)then ! precipitate! 334 d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0) 335 precip_rate(i) = precip_rate(i) - d_ql(i,k)*dmass(i,k)/ptimestep 336 endif 337 endif 338 enddo 339 340 elseif (precip_scheme_generic.ge.2) then 341 349 342 do i = 1, ngrid 350 343 if (rneb(i,k).GT.0.0) then 351 ! this is the ONLY place where zneb, precip_timescale_generic and cloud_sat_generic are used 352 353 zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale_generic)) * zoliq(i) & 354 * (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat_generic)**2)) * zfrac(i) 355 zrhol(i) = zrho(i) * zoliq(i) / zneb(i) 356 zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & 357 *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... 358 ztot(i) = zchau(i) + zfroi(i) 359 360 if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0 361 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) 362 zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) 363 endif 344 zoliq(i) = ql(i,k) 345 zrho(i) = pplay(i,k) / ( zt(i,k) * R ) 346 zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g) 347 zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds) 348 zfrac(i) = MAX(zfrac(i), 0.0) 349 zfrac(i) = MIN(zfrac(i), 1.0) 350 zneb(i) = MAX(rneb(i,k), seuil_neb) 351 endif 364 352 enddo 365 enddo 366 367 !precip scheme modified from Sundquist 78 (in q**3) 368 CASE(3) 369 do n = 1, ninter 353 354 !recalculate liquid GCS particle radii 355 call h2o_cloudrad(ngrid,nlayer,ql,reffh2oliq,reffh2oice) 356 357 SELECT CASE(precip_scheme_generic) 358 359 !precip scheme from Sundquist 78 360 CASE(2) 361 362 do n = 1, ninter 363 do i = 1, ngrid 364 if (rneb(i,k).GT.0.0) then 365 ! this is the ONLY place where zneb, precip_timescale_generic and cloud_sat_generic are used 366 367 zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale_generic)) * zoliq(i) & 368 * (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat_generic)**2)) * zfrac(i) 369 zrhol(i) = zrho(i) * zoliq(i) / zneb(i) 370 zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & 371 *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... 372 ztot(i) = zchau(i) + zfroi(i) 373 374 if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0 375 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) 376 zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) 377 endif 378 enddo 379 enddo 380 381 !precip scheme modified from Sundquist 78 (in q**3) 382 CASE(3) 383 do n = 1, ninter 384 do i = 1, ngrid 385 if (rneb(i,k).GT.0.0) then 386 ! this is the ONLY place where zneb, precip_timescale_generic and cloud_sat_generic are used 387 388 zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale_generic*cloud_sat_generic**2)) * (zoliq(i)/zneb(i))**3 389 zrhol(i) = zrho(i) * zoliq(i) / zneb(i) 390 zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & 391 *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... 392 ztot(i) = zchau(i) + zfroi(i) 393 394 if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0 395 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) 396 zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) 397 endif 398 enddo 399 enddo 400 401 !precip scheme modified from Boucher 95 402 CASE(4) 403 do n = 1, ninter 404 do i = 1, ngrid 405 if (rneb(i,k).GT.0.0) then 406 ! this is the ONLY place where zneb and c1 are used 407 zchau(i) = ptimestep/FLOAT(ninter) *c1* zrho(i) & 408 *(zoliq(i)/zneb(i))**2*reffh2oliq(i,k)*zneb(i)* zfrac(i) 409 zrhol(i) = zrho(i) * zoliq(i) / zneb(i) 410 zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & 411 *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... 412 ztot(i) = zchau(i) + zfroi(i) 413 414 if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0 415 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) 416 zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) 417 endif 418 enddo 419 enddo 420 421 END SELECT ! precip_scheme_generic 422 423 ! Change in cloud density and surface GCS values 370 424 do i = 1, ngrid 371 425 if (rneb(i,k).GT.0.0) then 372 ! this is the ONLY place where zneb, precip_timescale_generic and cloud_sat_generic are used 373 374 zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale_generic*cloud_sat_generic**2)) * (zoliq(i)/zneb(i))**3 375 zrhol(i) = zrho(i) * zoliq(i) / zneb(i) 376 zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & 377 *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... 378 ztot(i) = zchau(i) + zfroi(i) 379 380 if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0 381 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) 382 zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) 383 endif 384 enddo 385 enddo 386 387 !precip scheme modified from Boucher 95 388 CASE(4) 389 do n = 1, ninter 390 do i = 1, ngrid 391 if (rneb(i,k).GT.0.0) then 392 ! this is the ONLY place where zneb and c1 are used 393 zchau(i) = ptimestep/FLOAT(ninter) *c1* zrho(i) & 394 *(zoliq(i)/zneb(i))**2*reffh2oliq(i,k)*zneb(i)* zfrac(i) 395 zrhol(i) = zrho(i) * zoliq(i) / zneb(i) 396 zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & 397 *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... 398 ztot(i) = zchau(i) + zfroi(i) 399 400 if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0 401 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) 402 zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) 426 d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep 427 precip_rate(i) = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmass(i,k)/ptimestep 403 428 endif 404 429 enddo 405 enddo 406 407 END SELECT ! precip_scheme_generic 408 409 ! Change in cloud density and surface GCS values 430 431 endif ! if precip_scheme_generic=1 432 433 enddo ! of do k = nlayer, 1, -1 434 435 ! Rain or snow on the ground 410 436 do i = 1, ngrid 411 if (rneb(i,k).GT.0.0)then412 d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep413 precip_rate(i) = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmass(i,k)/ptimestep437 if(precip_rate(i).lt.0.0)then 438 print*,'Droplets of negative rain are falling...' 439 call abort 414 440 endif 441 if (t(i,1) .LT. T_h2O_ice_liq) then 442 dqssnow_generic(i,i_ice_generic) = precip_rate(i) 443 dqsrain_generic(i,i_ice_generic) = 0.0 444 else 445 dqssnow_generic(i,i_ice_generic) = 0.0 446 dqsrain_generic(i,i_ice_generic) = precip_rate(i) ! liquid water = ice for now 447 endif 448 449 ! For now we force : 450 dqsrain_generic(i,i_ice_generic) = precip_rate(i) 451 dqssnow_generic(i,i_ice_generic) = 0.0 415 452 enddo 416 453 417 endif ! if precip_scheme_generic=1 418 419 enddo ! of do k = nlayer, 1, -1 420 421 ! Rain or snow on the ground 422 do i = 1, ngrid 423 if(precip_rate(i).lt.0.0)then 424 print*,'Droplets of negative rain are falling...' 425 call abort 426 endif 427 if (t(i,1) .LT. T_h2O_ice_liq) then 428 dqssnow_generic(i) = precip_rate(i) 429 dqsrain_generic(i) = 0.0 430 else 431 dqssnow_generic(i) = 0.0 432 dqsrain_generic(i) = precip_rate(i) ! liquid water = ice for now 433 endif 434 enddo 435 436 ! now subroutine -----> GCM variables 437 if (evap_prec_generic) then 438 dqrain_generic(1:ngrid,1:nlayer,i_vap_generic)=d_q(1:ngrid,1:nlayer)/ptimestep 439 d_t(1:ngrid,1:nlayer)=d_t(1:ngrid,1:nlayer)/ptimestep 440 do i=1,ngrid 441 reevap_precip(i)=0. 442 do k=1,nlayer 443 reevap_precip(i)=reevap_precip(i)+dqrain_generic(i,k,i_vap_generic)*dmass(i,k) 444 enddo 445 enddo 446 else 447 dqrain_generic(1:ngrid,1:nlayer,i_vap_generic)=0.0 448 d_t(1:ngrid,1:nlayer)=0.0 449 endif 450 dqrain_generic(1:ngrid,1:nlayer,i_ice_generic) = d_ql(1:ngrid,1:nlayer)/ptimestep 454 ! now subroutine -----> GCM variables 455 if (evap_prec_generic) then 456 dq_rain_generic_vap(1:ngrid,1:nlayer,i_vap_generic)=d_q(1:ngrid,1:nlayer)/ptimestep 457 d_t(1:ngrid,1:nlayer)=d_t(1:ngrid,1:nlayer)/ptimestep 458 do i=1,ngrid 459 reevap_precip(i,i_vap_generic)=0. 460 do k=1,nlayer 461 reevap_precip(i,i_vap_generic)=reevap_precip(i,i_vap_generic)+dq_rain_generic_vap(i,k,i_vap_generic)*dmass(i,k) 462 enddo 463 enddo 464 else 465 dq_rain_generic_vap(1:ngrid,1:nlayer,i_vap_generic)=0.0 466 d_t(1:ngrid,1:nlayer)=0.0 467 endif 468 dq_rain_generic_cld(1:ngrid,1:nlayer,i_ice_generic) = d_ql(1:ngrid,1:nlayer)/ptimestep 469 470 end if ! if(one_call_ice_vap_generic) 471 472 end do ! do iq=1,nq loop on tracers 451 473 452 474 end subroutine rain_generic
Note: See TracChangeset
for help on using the changeset viewer.
