- Timestamp:
- Sep 1, 2020, 11:23:30 AM (4 years ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 5 edited
-
README (modified) (1 diff)
-
libf/phymars/callsedim_mod.F (modified) (2 diffs)
-
libf/phymars/growthrate.F (modified) (3 diffs)
-
libf/phymars/improvedclouds_mod.F (modified) (7 diffs)
-
libf/phymars/microphys.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r2404 r2407 3134 3134 a MCS variable among temp (default), dust and water ice, serving as reference for the binning. 3135 3135 More details in the preface of simu_MCS.F90 and in simu_MCS.def 3136 3137 == 01/09/2020 == MV 3138 Implementation of HDO in the microphysics of water ice clouds: 3139 - improvedclouds_mod.F: calculation of the HDO flux 3140 - growthrate.F, microphys.h: addings to take into account the effect of diffusion kinetics on fractionation 3141 - callsedim_mod.F: sedimentation of HDO as an isotope of water in the microphysics case -
trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F
r2332 r2407 516 516 & ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud, 517 517 & zqi(1,1,iq),wq,beta) 518 c Surface Tendencies519 c ~~~~~~~~~~520 do ig=1,ngrid521 pdqs_sed(ig,iq)=wq(ig,1)/ptimestep522 end do523 518 else 524 519 c water ice sedimentation … … 527 522 & ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq), 528 523 & zqi(1,1,iq),wq,beta) 529 c Surface tendencies 530 c ~~~~~~~~~~ 531 do ig=1,ngrid 532 pdqs_sed(ig,iq)=wq(ig,1)/ptimestep 524 endif ! of if (microphys) 525 c Surface tendencies 526 c ~~~~~~~~~~ 527 do ig=1,ngrid 528 pdqs_sed(ig,iq)=wq(ig,1)/ptimestep 529 end do 530 c Special case of isotopes 531 c ~~~~~~~~~~ 532 !MVals: for now only water can have an isotope, a son ("fils"), and it has to be hdo 533 if (nqfils(iq).gt.0) then 534 if (iq.eq.igcm_h2o_ice) then 535 iq2=igcm_hdo_ice 536 else 537 call abort_physic("callsedim_mod","invalid isotope",1) 538 endif 539 !MVals: input parameters in vlz_fi for hdo 540 do l=1,nlay 541 do ig=1,ngrid 542 if (zq0(ig,l,iq).gt.qperemin) then 543 Ratio0(ig,l)=zq0(ig,l,iq2)/zq0(ig,l,iq) 544 else 545 Ratio0(ig,l)=0. 546 endif 547 Ratio(ig,l)=Ratio0(ig,l) 548 masseq(ig,l)=max(masse(ig,l)*zq0(ig,l,iq),masseqmin) 549 w(ig,l)=wq(ig,l) !MVals: very important: hdo is transported by h2o (see vlsplt_p.F: correction bugg 15 mai 2015) 550 enddo !ig=1,ngrid 551 enddo !l=1,nlay 552 !MVals: no need to enter newsedim as the transporting mass w has been already calculated 553 call vlz_fi(ngrid,nlay,Ratio,2.,masseq,w,wq) 554 zqi(:,nlay,iq2)=zqi(:,nlay,iq)*Ratio0(:,nlay) 555 do l=1,nlay-1 556 do ig=1,ngrid 557 newmasse=max(masseq(ig,l)+w(ig,l+1)-w(ig,l),masseqmin) 558 Ratio(ig,l)=(Ratio0(ig,l)*masseq(ig,l) 559 & +wq(ig,l+1)-wq(ig,l))/newmasse 560 zqi(ig,l,iq2)=zqi(ig,l,iq)*Ratio(ig,l) 561 enddo 562 enddo !l=1,nlay-1 563 !MVals: hdo surface tendency 564 do ig=1,ngrid 565 if (w(ig,1).gt.masseqmin) then 566 pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*(wq(ig,1)/w(ig,1)) 567 else 568 pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*Ratio0(ig,1) 569 endif 533 570 end do 534 c Special case of isotopes 535 c ~~~~~~~~~~ 536 !MVals: Loop over the sons ("fils") 537 if (nqfils(iq).gt.0) then 538 if (iq.eq.igcm_h2o_ice) then 539 iq2=igcm_hdo_ice 540 else 541 call abort_physic("callsedim_mod","invalid isotope",1) 542 endif 543 !MVals: input paramters in vlz_fi for hdo 544 do l=1,nlay 545 do ig=1,ngrid 546 if (zq0(ig,l,iq).gt.qperemin) then 547 Ratio0(ig,l)=zq0(ig,l,iq2)/zq0(ig,l,iq) 548 else 549 Ratio0(ig,l)=0. 550 endif 551 Ratio(ig,l)=Ratio0(ig,l) 552 masseq(ig,l)=max(masse(ig,l)*zq0(ig,l,iq),masseqmin) 553 w(ig,l)=wq(ig,l) !MVals: very important: hdo is transported by h2o (see vlsplt_p.F: correction bugg 15 mai 2015) 554 enddo !ig=1,ngrid 555 enddo !l=1,nlay 556 !MVals: no need to enter newsedim as the transporting mass w has been already calculated 557 call vlz_fi(ngrid,nlay,Ratio,2.,masseq,w,wq) 558 zqi(:,nlay,iq2)=zqi(:,nlay,iq)*Ratio0(:,nlay) 559 do l=1,nlay-1 560 do ig=1,ngrid 561 newmasse=max(masseq(ig,l)+w(ig,l+1)-w(ig,l),masseqmin) 562 Ratio(ig,l)=(Ratio0(ig,l)*masseq(ig,l) 563 & +wq(ig,l+1)-wq(ig,l))/newmasse 564 zqi(ig,l,iq2)=zqi(ig,l,iq)*Ratio(ig,l) 565 enddo 566 enddo !l=1,nlay-1 567 !MVals: hdo surface tendency 568 do ig=1,ngrid 569 if (w(ig,1).gt.masseqmin) then 570 pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*(wq(ig,1)/w(ig,1)) 571 else 572 pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*Ratio0(ig,1) 573 endif 574 end do 575 endif !(nqfils(iq).gt.0) 576 endif ! of if (microphys) 577 571 endif !(nqfils(iq).gt.0) 578 572 c ----------------------------------------------------------------- 579 573 c GENERAL CASE -
trunk/LMDZ.MARS/libf/phymars/growthrate.F
r1266 r2407 1 subroutine growthrate(temp,pmid,psat,rcrystal,res )1 subroutine growthrate(temp,pmid,psat,rcrystal,res,Dv) 2 2 3 3 use tracer_mod, only: rho_ice … … 34 34 c Output 35 35 REAL res ! growth resistance (res=Rk+Rd) 36 36 REAL Dv ! water vapor diffusion coefficient 37 37 38 38 c local: … … 41 41 REAL k,Lv 42 42 REAL knudsen ! Knudsen number (gas mean free path/particle radius) 43 REAL afactor, Dv,lambda ! Intermediate computations for growth rate43 REAL afactor,lambda ! Intermediate computations for growth rate 44 44 REAL Rk,Rd 45 45 -
trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F
r2304 r2407 14 14 & igcm_h2o_ice, igcm_dust_mass, 15 15 & igcm_dust_number, igcm_ccn_mass, 16 & igcm_ccn_number 16 & igcm_ccn_number, 17 & igcm_hdo_vap,igcm_hdo_ice, 18 & qperemin 17 19 use conc_mod, only: mmean 18 20 use comcstfi_h, only: pi, cpp … … 82 84 REAL cste 83 85 REAL dMice ! mass of condensed ice 86 REAL dMice_hdo ! mass of condensed HDO ice 87 REAL alpha(ngrid,nlay) ! HDO equilibrium fractionation coefficient (Saturation=1) 88 REAL alpha_c(ngrid,nlay) ! HDO real fractionation coefficient 84 89 ! REAL sumcheck 85 90 REAL*8 ph2o ! Water vapor partial pressure (Pa) … … 103 108 104 109 REAL res ! Resistance growth 110 REAL Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient 105 111 106 112 … … 364 370 c get resistance growth 365 371 call growthrate(zt(ig,l),pplay(ig,l), 366 & real(ph2o/satu),rice(ig,l),res )372 & real(ph2o/satu),rice(ig,l),res,Dv) 367 373 368 374 res_out(ig,l) = res … … 391 397 subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep 392 398 393 399 c Special case of the isotope of water HDO 400 if (hdo) then 401 !! condensation 402 if (dMice.gt.0.0) then 403 !! do we use fractionation? 404 if (hdofrac) then 405 !! Calculation of the HDO vapor coefficient 406 Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) ) 407 & * kbz * zt(ig,l) / 408 & ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo) 409 & * sqrt(1.+mhdo/mco2) ) 410 !! Calculation of the fractionnation coefficient at equilibrium 411 alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2) 412 c alpha = exp(13525./zt(ig,l)**2.-5.59d-2) !Lamb 413 !! Calculation of the 'real' fractionnation coefficient 414 alpha_c(ig,l) = (alpha(ig,l)*satu)/ 415 & ( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.) 416 c alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics 417 else 418 alpha_c(ig,l) = 1.d0 419 endif 420 if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then 421 dMice_hdo= 422 & dMice*alpha_c(ig,l)* 423 & ( zq0(ig,l,igcm_hdo_vap) 424 & /zq0(ig,l,igcm_h2o_vap) ) 425 else 426 dMice_hdo=0. 427 endif 428 !! sublimation 429 else 430 if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then 431 dMice_hdo= 432 & dMice* 433 & ( zq0(ig,l,igcm_hdo_ice) 434 & /zq0(ig,l,igcm_h2o_ice) ) 435 else 436 dMice_hdo=0. 437 endif 438 endif !if (dMice.gt.0.0) 439 440 dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice)) 441 dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap)) 442 443 zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo 444 zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo 445 446 endif ! if (hdo) 394 447 395 448 !============================================================= … … 405 458 & + zq(ig,l,igcm_h2o_ice) 406 459 zq(ig,l,igcm_h2o_ice) = 0. 460 if (hdo) then 461 zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) 462 & + zq(ig,l,igcm_hdo_ice) 463 zq(ig,l,igcm_hdo_ice) = 0. 464 endif 407 465 c Dust particles 408 466 zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) … … 429 487 & (zq(1:ngrid,1:nlay,igcm_h2o_ice) - 430 488 & zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep 489 if (hdo) then 490 subpdqcloud(1:ngrid,1:nlay,igcm_hdo_vap) = 491 & (zq(1:ngrid,1:nlay,igcm_hdo_vap) - 492 & zq0(1:ngrid,1:nlay,igcm_hdo_vap))/microtimestep 493 subpdqcloud(1:ngrid,1:nlay,igcm_hdo_ice) = 494 & (zq(1:ngrid,1:nlay,igcm_hdo_ice) - 495 & zq0(1:ngrid,1:nlay,igcm_hdo_ice))/microtimestep 496 endif 431 497 subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) = 432 498 & (zq(1:ngrid,1:nlay,igcm_ccn_mass) - -
trunk/LMDZ.MARS/libf/phymars/microphys.h
r1816 r2407 18 18 ! Molecular weight of H2O (kg.mol-1) 19 19 DOUBLE PRECISION, PARAMETER :: mh2o = 18.d-3 20 ! Molecular weight of HDO (kg.mol-1) 21 DOUBLE PRECISION, PARAMETER :: mhdo = 19.d-3 20 22 ! Molecular weight of CO2 (kg.mol-1) 21 23 DOUBLE PRECISION, PARAMETER :: mco2 = 44.d-3 … … 27 29 ! Effective H2O gas molecular radius (m) 28 30 DOUBLE PRECISION, PARAMETER :: molh2o = 1.2d-10 31 ! Effective HDO gas molecular radius (m) 32 DOUBLE PRECISION, PARAMETER :: molhdo = 1.2d-10 29 33 ! Surface tension of ice/vapor (N.m) 30 34 DOUBLE PRECISION, PARAMETER :: sigh2o = 0.12
Note: See TracChangeset
for help on using the changeset viewer.
