Changeset 1685
- Timestamp:
- Apr 3, 2017, 11:08:39 AM (8 years ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 8 edited
-
README (modified) (1 diff)
-
libf/phymars/co2cloud.F (modified) (32 diffs)
-
libf/phymars/improvedCO2clouds.F (modified) (33 diffs)
-
libf/phymars/initracer.F (modified) (2 diffs)
-
libf/phymars/massflowrateCO2.F (modified) (1 diff)
-
libf/phymars/microphys.h (modified) (4 diffs)
-
libf/phymars/nucleaCO2.F (modified) (6 diffs)
-
libf/phymars/physiq_mod.F (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r1684 r1685 2427 2427 solvarmod==0. 2428 2428 Guidelines for min/ave/max EUV input: fixed_euv_value=80/140/320. 2429 2430 2431 == 03/04/2017 == JA 2432 further work on the CO2 clouds scheme: 2433 -water ice clouds can now serve as condensation nucleii for CO2 clouds. Memory of CO2 CCN origin is kept, but a study of the water cycle is needed to know ig further work is necessary. 2434 -co2cloud.F, improvedCO2clouds.F, nucleaCO2.F, massflowrateCO2.F and physiq_mod.F have been modified accordingly. 2435 -co2 cloud scheme tuning has been improved for more stability. 2436 -a new contact parameter for CO2 ice on silicate is used (mtetaco2=0.78) in microphys.h. Reference in microphysic.h 2437 -CO2 ice temperature is computed as a function of temperature in CO2cloud.F and improvedCO2cloud.F. Reference is in initracer.F -
trunk/LMDZ.MARS/libf/phymars/co2cloud.F
r1655 r1685 1 1 SUBROUTINE co2cloud(ngrid,nlay,ptimestep, 2 2 & pplev,pplay,pdpsrf,pzlay,pt,pdt, 3 3 & pq,pdq,pdqcloudco2,pdtcloudco2, 4 4 & nq,tau,tauscaling,rdust,rice,riceco2,nuice, 5 & rsedcloudco2,rhocloudco2,zlev,pdqs_sedco2) 5 & rsedcloudco2,rhocloudco2, 6 & rsedcloud,rhocloud,zlev,pdqs_sedco2) 6 7 ! to use 'getin' 7 8 use dimradmars_mod, only: naerkind … … 11 12 use conc_mod, only: mmean 12 13 use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice, 13 & igcm_dust_mass, igcm_dust_number, 14 & igcm_dust_mass, igcm_dust_number,igcm_h2o_ice, 15 & igcm_ccn_mass,igcm_ccn_number, 14 16 & igcm_ccnco2_mass, igcm_ccnco2_number, 15 17 & rho_dust, nuiceco2_sed, nuiceco2_ref, 16 & rho_ice_co2,r3n_q, 18 & rho_ice_co2,r3n_q,rho_ice,nuice_sed, 17 19 & meteo_flux_mass,meteo_flux_number, 18 20 & meteo_alt … … 73 75 real rice(ngrid,nlay) ! Water Ice mass mean radius (m) 74 76 ! used for nucleation of CO2 on ice-coated ccns 75 77 DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) 76 78 REAL tau(ngrid,naerkind) ! Column dust optical depth at each point 77 79 REAL tauscaling(ngrid) ! Convertion factor for dust amount … … 95 97 c local: 96 98 c ------ 97 99 !water 100 real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius 101 real rhocloud(ngrid,nlay) ! Cloud density (kg.m-3) 98 102 ! for ice radius computation 99 103 REAL Mo,No … … 127 131 real epaisseur (ngrid,nlay) ! Layer thickness (m) 128 132 real masse (ngrid,nlay) ! Layer mass (kg.m-2) 129 130 133 double precision diff,diff0 131 134 integer meteo_lvl … … 134 137 real sav_trac(ngrid,nlay,nq) 135 138 real pdqsed(ngrid,nlay,nq) 139 140 DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) 141 DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:) 142 DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:) 136 143 c ** un petit test de coherence 137 144 c -------------------------- … … 144 151 stop 145 152 endif 146 153 write(*,*) "co2cloud.F: rho_ice_co2 = ",rho_ice_co2 154 147 155 write(*,*) "co2cloud: igcm_co2=",igcm_co2 148 156 write(*,*) " igcm_co2_ice=",igcm_co2_ice … … 161 169 write(*,*)"CO2 Microphysics timestep is",microtimestep 162 170 163 171 172 if (meteo_flux_number .ne. 0) then 173 164 174 write(*,*) "Warning ! Meteoritic flux of dust is turned on" 165 175 write(*,*) "Dust mass = ",meteo_flux_mass … … 167 177 write(*,*) "Are added at the z-level (km)",meteo_alt 168 178 write(*,*) "Every timestep in co2cloud.F" 169 179 endif 180 if (.not. allocated(memdMMccn)) allocate(memdMMccn(ngrid,nlay)) 181 if (.not. allocated(memdNNccn)) allocate(memdNNccn(ngrid,nlay)) 182 if (.not. allocated(memdMMh2o)) allocate(memdMMh2o(ngrid,nlay)) 183 184 memdMMccn(:,:)=0. 185 memdMMh2o(:,:)=0. 186 memdNNccn(:,:)=0. 170 187 firstcall=.false. 171 188 ENDIF ! of IF (firstcall) … … 173 190 c-----Initialization 174 191 192 beta=0.85 193 subpdq(1:ngrid,1:nlay,1:nq) = 0 194 subpdt(1:ngrid,1:nlay) = 0 195 subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0 196 subpdtcloudco2(1:ngrid,1:nlay) = 0 197 175 198 c add meteoritic flux of dust 176 199 !convert meteo_alt (in km) to z-level 177 200 !pzlay altitudes of the layers 178 201 if (meteo_flux_number .ne. 0) then 179 202 do ig=1, ngrid 180 203 diff0=1000 … … 188 211 enddo 189 212 c write(*,*) "meteo_lvl=",meteo_lvl 190 pq(ig,meteo_lvl,igcm_dust_mass)=pq(ig,meteo_lvl,igcm_dust_mass) 191 & +dble(meteo_flux_mass) 192 pq(ig,meteo_lvl,igcm_dust_number)= 193 & pq(ig,meteo_lvl,igcm_dust_number) 194 & +dble(meteo_flux_number) 213 pdq(ig,meteo_lvl,igcm_dust_mass)= 214 & pdq(ig,meteo_lvl,igcm_dust_mass) 215 & +dble(meteo_flux_mass)/tauscaling(ig) 216 pdq(ig,meteo_lvl,igcm_dust_number)= 217 & pdq(ig,meteo_lvl,igcm_dust_number) 218 & +dble(meteo_flux_number)/tauscaling(ig) 195 219 enddo 196 197 call WRITEDIAGFI(ngrid,"pzlay","altitude","km",3, 198 & pzlay) 199 beta=0.85 200 subpdq(1:ngrid,1:nlay,1:nq) = 0 201 subpdt(1:ngrid,1:nlay) = 0 202 subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0 203 subpdtcloudco2(1:ngrid,1:nlay) = 0 220 endif 221 c call WRITEDIAGFI(ngrid,"pzlay","altitude","km",3, 222 c & pzlay) 204 223 205 224 … … 220 239 masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 221 240 epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l) 222 223 241 enddo 224 242 enddo 225 226 227 228 229 230 243 231 244 … … 274 287 & subpdq(ig,l,igcm_co2_ice) 275 288 & + pdq(ig,l,igcm_co2_ice) 289 276 290 subpdq(ig,l,igcm_co2) = 277 291 & subpdq(ig,l,igcm_co2) 278 292 & + pdq(ig,l,igcm_co2) 293 294 subpdq(ig,l,igcm_h2o_ice) = 295 & subpdq(ig,l,igcm_h2o_ice) 296 & + pdq(ig,l,igcm_h2o_ice) 297 298 subpdq(ig,l,igcm_ccn_mass) = 299 & subpdq(ig,l,igcm_ccn_mass) 300 & + pdq(ig,l,igcm_ccn_mass) 301 302 subpdq(ig,l,igcm_ccn_number) = 303 & subpdq(ig,l,igcm_ccn_number) 304 & + pdq(ig,l,igcm_ccn_number) 279 305 280 306 tempo_traceur_t(ig,l)= pt(ig,l)+subpdt(ig,l)*microtimestep … … 289 315 DO l=1, nlay 290 316 DO ig=1,ngrid 317 318 rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* 319 & tempo_traceur_t(ig,l)-2.87e-6* 320 & tempo_traceur_t(ig,l)*tempo_traceur_t(ig,l)) 321 322 rho_ice_co2=rho_ice_co2T(ig,l) 291 323 Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30) 292 324 Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number), … … 296 328 mdustJA= tempo_traceurs(ig,l,igcm_dust_mass) 297 329 ndustJA=tempo_traceurs(ig,l,igcm_dust_number) 298 if (( ndustJA .lt. tauscaling(ig)) .or. (mdustJA.lt.330 if ((Nccnco2 .lt. tauscaling(ig)) .or. (Qccnco2 .lt. 299 331 & 1.e-30 *tauscaling(ig))) then 300 332 rdust(ig,l)=1.e-10 301 333 else 302 rdust(ig,l)=(3./4./pi/2500.* mdustJA/ndustJA)**(1./3.)303 rdust(ig,l)=m in(rdust(ig,l),5.e-6)304 rdust(ig,l)=max(rdust(ig,l),1.e-9)334 rdust(ig,l)=(3./4./pi/2500.*Qccnco2/Nccnco2)**(1./3.) 335 rdust(ig,l)=max(rdust(ig,l),1.e-10) 336 ! rdust(ig,l)=min(rdust(ig,l),5.e-6) 305 337 end if 306 crhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2307 c & + Qccnco2*rho_dust)308 c & / (Niceco2 + Qccnco2)309 riceco2(ig,l)= Niceco2*3.0/310 & (4.0*rho_ice_co2*pi*Nccnco2)311 & +rdust(ig,l)*rdust(ig,l)*rdust(ig,l)312 riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)338 rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 339 & + Qccnco2*tauscaling(ig)*rho_dust) 340 & / (Niceco2 + Qccnco2*tauscaling(ig)) 341 c riceco2(ig,l)= Niceco2*3.0/ 342 c & (4.0*rho_ice_co2*pi*Nccnco2) 343 c & +rdust(ig,l)*rdust(ig,l)*rdust(ig,l) 344 c riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0) 313 345 c write(*,*) "in co2clouds, rice = ",riceco2(ig,l) 314 346 c write(*,*) "in co2clouds, rho = ",rhocloudco2t(ig,l) 315 347 316 call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, 317 & tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l)) 348 riceco2(ig,l)=(Niceco2*3.0/ 349 & (4.0*rho_ice_co2*pi*Nccnco2 350 & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) 351 & *rdust(ig,l))**(1.0/3.0) 352 353 riceco2(ig,l)=max(1.e-10,riceco2(ig,l)) 354 !riceco2(ig,l)=min(5.e-5,riceco2(ig,l)) 355 356 c call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, 357 c & tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l)) 318 358 c write(*,*) "in co2clouds, rice update = ",riceco2(ig,l) 319 359 c write(*,*) "in co2clouds, rho update = " … … 323 363 & (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed), 324 364 & rdust(ig,l)) 325 rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l), 5.e-4)365 rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5) 326 366 c write(*,*) 'Rsedcloud = ',rsedcloudco2(ig,l) 327 367 !write(*,*) 'Rhocloudco2 = ',rhocloudco2t(ig,l) … … 411 451 & pplay,pt,subpdt, 412 452 & pq,subpdq,subpdqcloudco2,subpdtcloudco2, 413 & nq,tauscaling )453 & nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn) 414 454 415 455 ELSE 416 456 417 457 write(*,*) ' no simpleCO2clouds procedure: STOP' ! listo 418 STOP 419 420 c CALL simpleclouds(ngrid,nlay,microtimestep, ! for water-ice clouds 421 c & pplay,pzlay,pt,subpdt, 422 c & pq,subpdq,subpdqcloud,subpdtcloud, 423 c & nq,tau,riceco2) 458 write(*,*) ' microphysco2 and co2clouds must be .true.' ! listo 459 STOP 460 424 461 ENDIF 425 462 … … 458 495 & subpdq(ig,l,igcm_co2) 459 496 & + subpdqcloudco2(ig,l,igcm_co2) 497 498 subpdq(ig,l,igcm_h2o_ice) = 499 & subpdq(ig,l,igcm_h2o_ice) 500 & + subpdqcloudco2(ig,l,igcm_h2o_ice) 501 502 subpdq(ig,l,igcm_ccn_mass) = 503 & subpdq(ig,l,igcm_ccn_mass) 504 & + subpdqcloudco2(ig,l,igcm_ccn_mass) 505 506 subpdq(ig,l,igcm_ccn_number) = 507 & subpdq(ig,l,igcm_ccn_number) 508 & + subpdqcloudco2(ig,l,igcm_ccn_number) 460 509 ENDDO 461 510 ENDDO 462 511 463 464 !ici 465 ! call WRITEdiagfi(ngrid,"co2cloud000","co2 traceur","kg/kg",1, 466 ! & pq(1,:,igcm_co2_ice) + ptimestep* 467 ! & ( subpdq(1,:,igcm_co2_ice))) 512 468 513 469 514 … … 492 537 DO ig=1,ngrid 493 538 pdtcloudco2(ig,l) = 494 & subpdt(ig,l)/ imicro-pdt(ig,l)539 & subpdt(ig,l)/real(imicro)-pdt(ig,l) 495 540 ENDDO 496 541 ENDDO 497 542 498 543 c------ Tracers tendencies pdqcloud 499 544 DO l=1,nlay … … 501 546 502 547 pdqcloudco2(ig,l,igcm_co2_ice) = 503 & subpdq(ig,l,igcm_co2_ice)/ imicro548 & subpdq(ig,l,igcm_co2_ice)/real(imicro) 504 549 & - pdq(ig,l,igcm_co2_ice) 550 505 551 pdqcloudco2(ig,l,igcm_co2) = 506 & subpdq(ig,l,igcm_co2)/ imicro552 & subpdq(ig,l,igcm_co2)/real(imicro) 507 553 & - pdq(ig,l,igcm_co2) 554 555 pdqcloudco2(ig,l,igcm_h2o_ice) = 556 & subpdq(ig,l,igcm_h2o_ice)/real(imicro) 557 & - pdq(ig,l,igcm_h2o_ice) 508 558 ENDDO 509 559 ENDDO … … 515 565 516 566 517 IF(microphysco2) THEN567 c IF(microphysco2) THEN 518 568 DO l=1,nlay 519 569 DO ig=1,ngrid 520 570 pdqcloudco2(ig,l,igcm_ccnco2_mass) = 521 & subpdq(ig,l,igcm_ccnco2_mass)/ imicro571 & subpdq(ig,l,igcm_ccnco2_mass)/real(imicro) 522 572 & - pdq(ig,l,igcm_ccnco2_mass) 523 573 pdqcloudco2(ig,l,igcm_ccnco2_number) = 524 & subpdq(ig,l,igcm_ccnco2_number)/ imicro574 & subpdq(ig,l,igcm_ccnco2_number)/real(imicro) 525 575 & - pdq(ig,l,igcm_ccnco2_number) 576 pdqcloudco2(ig,l,igcm_ccn_mass) = 577 & subpdq(ig,l,igcm_ccn_mass)/real(imicro) 578 & - pdq(ig,l,igcm_ccn_mass) 579 pdqcloudco2(ig,l,igcm_ccn_number) = 580 & subpdq(ig,l,igcm_ccn_number)/real(imicro) 581 & - pdq(ig,l,igcm_ccn_number) 526 582 ENDDO 527 583 ENDDO 528 ENDIF584 c ENDIF 529 585 530 586 … … 553 609 & ptimestep* (pdq(ig,l,igcm_ccnco2_number) + 554 610 & pdqcloudco2(ig,l,igcm_ccnco2_number)) 555 & .lt. 0)556 & .or. (pq(ig,l,igcm_ccnco2_mass) +611 & .lt. 1.e-15) 612 & .or. (pq(ig,l,igcm_ccnco2_mass) + 557 613 & ptimestep* (pdq(ig,l,igcm_ccnco2_mass) + 558 614 & pdqcloudco2(ig,l,igcm_ccnco2_mass)) 559 & .lt. 0)) THEN615 & .lt. 1.e-30)) THEN 560 616 561 617 pdqcloudco2(ig,l,igcm_ccnco2_number) = 562 618 & - pq(ig,l,igcm_ccnco2_number)/ptimestep 563 & - pdq(ig,l,igcm_ccnco2_number) +0619 & - pdq(ig,l,igcm_ccnco2_number)+1.e-15 564 620 565 621 pdqcloudco2(ig,l,igcm_dust_number) = … … 568 624 pdqcloudco2(ig,l,igcm_ccnco2_mass) = 569 625 & - pq(ig,l,igcm_ccnco2_mass)/ptimestep 570 & - pdq(ig,l,igcm_ccnco2_mass)+ 0626 & - pdq(ig,l,igcm_ccnco2_mass)+1.e-30 571 627 572 628 pdqcloudco2(ig,l,igcm_dust_mass) = … … 584 640 IF ( (pq(ig,l,igcm_dust_number) + 585 641 & ptimestep* (pdq(ig,l,igcm_dust_number) + 586 & pdqcloudco2(ig,l,igcm_dust_number)) .l t. 0.)642 & pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.e-15) 587 643 & .or. (pq(ig,l,igcm_dust_mass)+ 588 644 & ptimestep* (pdq(ig,l,igcm_dust_mass) + 589 645 & pdqcloudco2(ig,l,igcm_dust_mass)) 590 & .l t. 0.)) then646 & .le. 1.e-30)) then 591 647 592 648 pdqcloudco2(ig,l,igcm_dust_number) = 593 649 & - pq(ig,l,igcm_dust_number)/ptimestep 594 & - pdq(ig,l,igcm_dust_number)+ 0650 & - pdq(ig,l,igcm_dust_number)+1.e-15 595 651 596 652 pdqcloudco2(ig,l,igcm_ccnco2_number) = … … 599 655 pdqcloudco2(ig,l,igcm_dust_mass) = 600 656 & - pq(ig,l,igcm_dust_mass)/ptimestep 601 & - pdq(ig,l,igcm_dust_mass) + 0657 & - pdq(ig,l,igcm_dust_mass) +1.e-30 602 658 603 659 pdqcloudco2(ig,l,igcm_ccnco2_mass) = … … 609 665 610 666 611 DO l=1,nlay 612 DO ig=1,ngrid 613 IF (pq(ig,l,igcm_co2_ice) + ptimestep* 614 & (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) 615 & .lt. 1.e-25) THEN 616 pdqcloudco2(ig,l,igcm_co2_ice) = 617 & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice) 618 & +1.e-25 619 pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice) 620 ENDIF 621 ENDDO 622 ENDDO 667 c DO l=1,nlay 668 c DO ig=1,ngrid 669 c IF (pq(ig,l,igcm_co2_ice) + ptimestep* 670 c & (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) 671 c & .lt. 1.e-25) THEN 672 c pdqcloudco2(ig,l,igcm_co2_ice) = 673 c & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice) 674 c pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice) 675 c write(*,*) "WARNING CO2 ICE in co2cloud.F" 676 677 c ENDIF 678 c IF (pq(ig,l,igcm_co2) + ptimestep* 679 c & (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2)) 680 c & .lt. 1.e-10) THEN 681 c pdqcloudco2(ig,l,igcm_co2) = 682 c & - pq(ig,l,igcm_co2)/ptimestep - pdq(ig,l,igcm_co2) 683 c pdqcloudco2(ig,l,igcm_co2_ice) = -pdqcloudco2(ig,l,igcm_co2) 684 c write(*,*) "WARNING CO2 in co2cloud.F" 685 c ENDIF 686 c ENDDO 687 c ENDDO 623 688 624 689 … … 628 693 c------Only rsedcloudco2 is used for the co2 (cloud) cycle 629 694 630 IF(scavenging) THEN631 DO l=1, nlay632 DO ig=1,ngrid695 c IF(scavenging) THEN 696 c DO l=1, nlay 697 c DO ig=1,ngrid 633 698 634 699 c call updaterdust( … … 641 706 c & rdust(ig,l)) 642 707 c write(*,*) "in co2clouds, rdust(ig,l)= ",rdust(ig,l) 643 mdustJA= pq(ig,l,igcm_dust_mass) +644 & (pdq(ig,l,igcm_dust_mass) +645 & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep646 ndustJA=pq(ig,l,igcm_dust_number) +647 & (pdq(ig,l,igcm_dust_number) +648 & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep649 if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt.650 & 1.e-30 *tauscaling(ig))) then651 rdust(ig,l)=1.e-10652 else653 rdust(ig,l)=(3./4./pi/2500.*mdustJA/ndustJA)**(1./3.)654 rdust(ig,l)=min(rdust(ig,l),5.e-4)655 rdust(ig,l)=max(rdust(ig,l),1.e-10)656 endif657 ENDDO658 ENDDO659 ENDIF708 c mdustJA= pq(ig,l,igcm_dust_mass) + 709 c & (pdq(ig,l,igcm_dust_mass) + 710 c & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep 711 c ndustJA=pq(ig,l,igcm_dust_number) + 712 c & (pdq(ig,l,igcm_dust_number) + 713 c & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep 714 c if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt. 715 c & 1.e-30 *tauscaling(ig))) then 716 c rdust(ig,l)=1.e-10 717 c else 718 c rdust(ig,l)=(3./4./pi/2500.*mdustJA/ndustJA)**(1./3.) 719 c rdust(ig,l)=min(rdust(ig,l),5.e-6) 720 c rdust(ig,l)=max(rdust(ig,l),1.e-10) 721 c endif 722 c ENDDO 723 c ENDDO 724 c ENDIF 660 725 661 726 … … 690 755 & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)* 691 756 & tauscaling(ig),1.e-30) 692 c rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 + Qccnco2*rho_dust) 693 c & / (Niceco2 + Qccnco2) 694 c rhocloudco2(ig,l) = min(max(rhocloudco2t,rho_ice_co2),rho_dust) 695 696 c write(*,*) "test, nccnco2 =",nccnco22 697 698 699 riceco2(ig,l)= Niceco2*3.0/ 757 758 759 if ((Nccnco2 .lt. 1) .or. (Qccnco2 .lt. 760 & 1.e-30)) then 761 rdust(ig,l)=1.e-10 762 else 763 764 rdust(ig,l)= Qccnco2 765 & *0.75/pi/rho_dust 766 & / Nccnco2 767 rdust(ig,l)= rdust(ig,l)**(1./3.) 768 rdust(ig,l)=max(1.e-10,rdust(ig,l)) 769 !rdust(ig,l)=min(5.e-5,rdust(ig,l)) 770 endif 771 772 rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* 773 & (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep) 774 & -2.87e-6* 775 & (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep)* 776 & (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep)) 777 rho_ice_co2=rho_ice_co2T(ig,l) 778 779 riceco2(ig,l)=(Niceco2*3.0/ 700 780 & (4.0*rho_ice_co2*pi*Nccnco2) 701 & +rdust(ig,l)*rdust(ig,l)*rdust(ig,l) 702 703 riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0) 704 c write(*,*) "In co2cloud, after loop, riceco2 =",riceco2(ig,l) 705 c write(*,*) "In co2cloud, after loop, rhoco2 =" 706 c & ,rhocloudco2t(ig,l) 707 708 call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, 709 & tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l)) 710 711 c write(*,*) "In co2cloud, after loop and update, riceco2 =" 712 c & ,riceco2(ig,l) 713 c write(*,*) "In co2cloud, after loop and update, rhoco2 =" 714 c & ,rhocloudco2t(ig,l) 715 716 if ( Niceco2 717 & .le. 1.e-23 .or. riceco2(ig,l) .le. 1.e-10 .or. 718 & riceco2(ig,l) .ge. 4.999e-4) then ! .or. riceco2(ig,l) .gt. 1.e-4 ) then 781 & +rdust(ig,l)*rdust(ig,l) 782 & *rdust(ig,l) )**(1.0/3.0) 783 784 if ( (Niceco2 .le. 1.e-23) .or. 785 & (riceco2(ig,l) .le. rdust(ig,l)) 786 & .or. (Nccnco2 .le. 0.01))then 787 788 c write(*,*) "Riceco2 co2cloud before reset=",riceco2(ig,l) 789 c write(*,*) "Niceco2 co2cloud before reset=",Niceco2 790 c write(*,*) "Nccnco2 co2cloud before reset=",Nccnco2 791 c write(*,*) "Rdust co2cloud before reset=",rdust(ig,l) 792 793 c$$$ !NO CLOUD : RESET TRACERS AND CONSERVE MASS 794 if (pq(ig,l,igcm_co2_ice)+(pdq(ig,l,igcm_co2_ice)+ 795 & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep .le. 0.) then 796 pdqcloudco2(ig,l,igcm_co2)=0. 797 pdqcloudco2(ig,l,igcm_co2_ice)=0. 798 else 799 pdqcloudco2(ig,l,igcm_co2)= pq(ig,l,igcm_co2_ice) 800 & /ptimestep+pdq(ig,l,igcm_co2_ice) 801 pdqcloudco2(ig,l,igcm_co2_ice)=-pq(ig,l,igcm_co2_ice) 802 & /ptimestep-pdq(ig,l,igcm_co2_ice) 803 endif 804 !Reverse h2o ccn and ices into h2o tracers 805 if (memdMMccn(ig,l) .gt. 0) then 806 pdqcloudco2(ig,l,igcm_ccn_mass)=memdMMccn(ig,l) 807 else 808 memdMMccn(ig,l)=0. 809 pdqcloudco2(ig,l,igcm_ccn_mass)=0. 810 endif 811 if (memdNNccn(ig,l) .gt. 0) then 812 pdqcloudco2(ig,l,igcm_ccn_number)=memdNNccn(ig,l) 813 else 814 memdNNccn(ig,l)=0. 815 pdqcloudco2(ig,l,igcm_ccn_number)=0. 816 endif 817 if (memdMMh2o(ig,l) .gt. 0) then 818 pdqcloudco2(ig,l,igcm_h2o_ice)=memdMMh2o(ig,l) 819 else 820 memdMMh2o(ig,l)=0. 821 pdqcloudco2(ig,l,igcm_h2o_ice)=0. 822 endif 823 824 if (pq(ig,l,igcm_ccnco2_mass)+(pdq(ig,l,igcm_ccnco2_mass)+ 825 & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep 826 & .le. 1.e-30) then 827 pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. 828 pdqcloudco2(ig,l,igcm_ccnco2_number)=0. 829 pdqcloudco2(ig,l,igcm_co2)=0. 830 pdqcloudco2(ig,l,igcm_co2_ice)=0. 831 else 832 pdqcloudco2(ig,l,igcm_ccnco2_mass)= 833 & -pq(ig,l,igcm_ccnco2_mass) 834 & /ptimestep-pdq(ig,l,igcm_ccnco2_mass) 835 endif 836 837 if (pq(ig,l,igcm_ccnco2_number)+ 838 & (pdq(ig,l,igcm_ccnco2_number)+ 839 & pdqcloudco2(ig,l,igcm_ccnco2_number)) 840 & *ptimestep.le. 1.e-30) 841 & then 842 pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. 843 pdqcloudco2(ig,l,igcm_ccnco2_number)=0. 844 pdqcloudco2(ig,l,igcm_co2)=0. 845 pdqcloudco2(ig,l,igcm_co2_ice)=0. 846 else 847 pdqcloudco2(ig,l,igcm_ccnco2_number)= 848 & -pq(ig,l,igcm_ccnco2_number) 849 & /ptimestep-pdq(ig,l,igcm_ccnco2_number) 850 endif 851 852 if (pq(ig,l,igcm_dust_number)+ 853 & (pdq(ig,l,igcm_dust_number)+ 854 & pdqcloudco2(ig,l,igcm_dust_number)) 855 & *ptimestep.le. 1.e-30) 856 & then 857 pdqcloudco2(ig,l,igcm_dust_number)=0. 858 pdqcloudco2(ig,l,igcm_dust_mass)=0. 859 else 860 pdqcloudco2(ig,l,igcm_dust_number)= 861 & pq(ig,l,igcm_ccnco2_number) 862 & /ptimestep+pdq(ig,l,igcm_ccnco2_number) 863 & -memdNNccn(ig,l) 864 endif 865 866 if (pq(ig,l,igcm_dust_mass)+ 867 & (pdq(ig,l,igcm_dust_mass)+ 868 & pdqcloudco2(ig,l,igcm_dust_mass)) 869 & *ptimestep .le. 1.e-30) 870 & then 871 pdqcloudco2(ig,l,igcm_dust_number)=0. 872 pdqcloudco2(ig,l,igcm_dust_mass)=0. 873 else 874 pdqcloudco2(ig,l,igcm_dust_mass)= 875 & pq(ig,l,igcm_ccnco2_mass) 876 & /ptimestep+pdq(ig,l,igcm_ccnco2_mass) 877 & -(memdMMh2o(ig,l)+memdMMccn(ig,l)) 878 endif 879 memdMMccn(ig,l)=0. 880 memdMMh2o(ig,l)=0. 881 memdNNccn(ig,l)=0. 719 882 riceco2(ig,l)=0. 883 pdtcloudco2(ig,l)=0. 884 endif 885 886 720 887 721 !NO CLOUD : RESET TRACER AND CONSERVE MASS 722 pdqcloudco2(ig,l,igcm_co2)= pq(ig,l,igcm_co2_ice) 723 & /ptimestep+pdq(ig,l,igcm_co2_ice) 724 725 pdqcloudco2(ig,l,igcm_co2_ice)=-pq(ig,l,igcm_co2_ice) 726 & /ptimestep-pdq(ig,l,igcm_co2_ice) 727 728 pdqcloudco2(ig,l,igcm_ccnco2_mass)= 729 & -pq(ig,l,igcm_ccnco2_mass) 730 & /ptimestep-pdq(ig,l,igcm_ccnco2_mass) 731 732 pdqcloudco2(ig,l,igcm_ccnco2_number)= 733 & -pq(ig,l,igcm_ccnco2_number) 734 & /ptimestep-pdq(ig,l,igcm_ccnco2_number) 735 736 pdqcloudco2(ig,l,igcm_dust_number)= 737 & pq(ig,l,igcm_ccnco2_number) 738 & /ptimestep+pdq(ig,l,igcm_ccnco2_number) 739 740 pdqcloudco2(ig,l,igcm_dust_mass)= 741 & pq(ig,l,igcm_ccnco2_mass) 742 & /ptimestep+pdq(ig,l,igcm_ccnco2_mass) 743 c$$$ 744 745 c$$$ 746 endif 747 748 c write(*,*) "in co2clouds, riceco2(ig,l)v2= ",riceco2(ig,l) 749 750 ENDDO 751 ENDDO 752 888 !update rice water 889 call updaterice_micro( 890 & pq(ig,l,igcm_h2o_ice) + ! ice mass 891 & (pdq(ig,l,igcm_h2o_ice) + ! ice mass 892 & pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass 893 & pq(ig,l,igcm_ccn_mass) + ! ccn mass 894 & (pdq(ig,l,igcm_ccn_mass) + ! ccn mass 895 & pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, ! ccn mass 896 & pq(ig,l,igcm_ccn_number) + ! ccn number 897 & (pdq(ig,l,igcm_ccn_number) + ! ccn number 898 & pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number 899 & tauscaling(ig),rice(ig,l),rhocloud(ig,l)) 900 901 902 call updaterdust( 903 & pq(ig,l,igcm_dust_mass) + ! dust mass 904 & (pdq(ig,l,igcm_dust_mass) + ! dust mass 905 & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass 906 & pq(ig,l,igcm_dust_number) + ! dust number 907 & (pdq(ig,l,igcm_dust_number) + ! dust number 908 & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number 909 & rdust(ig,l)) 910 911 ENDDO 912 ENDDO 913 914 915 753 916 ELSE ! no microphys ! not of concern for co2 clouds - listo 754 917 … … 758 921 c TO CHEK for relevancy - listo 759 922 923 c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005 924 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 925 c Then that should not affect the ice particle radius 926 !do ig=1,ngrid 927 ! if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then 928 ! if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3))) 929 &! riceco2(ig,2)=riceco2(ig,3) 930 ! riceco2(ig,1)=riceco2(ig,2) 931 ! end if 932 !end do 933 760 934 c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005 761 935 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ … … 772 946 DO l=1,nlay 773 947 DO ig=1,ngrid 948 rsedcloud(ig,l)=max(rice(ig,l)* 949 & (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed), 950 & rdust(ig,l)) 951 ! rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4) 952 ENDDO 953 ENDDO 954 955 DO l=1,nlay 956 DO ig=1,ngrid 774 957 rsedcloudco2(ig,l)=max(riceco2(ig,l)* 775 958 & (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed), 776 959 & rdust(ig,l)) 777 rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e- 4)960 rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5) 778 961 ENDDO 779 962 ENDDO … … 782 965 do ig=1,ngrid 783 966 do l=1,nlay 784 satuco2(ig,l) = pq(ig,l,igcm_co2)* 967 satuco2(ig,l) = (pq(ig,l,igcm_co2) + 968 & (pdq(ig,l,igcm_co2) + 969 & pdqcloudco2(ig,l,igcm_co2))*ptimestep)* 785 970 & (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l) 786 971 -
trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F
r1649 r1685 2 2 & pplay,pt,pdt, 3 3 & pq,pdq,pdqcloudco2,pdtcloudco2, 4 & nq,tauscaling) 4 & nq,tauscaling, 5 & memdMMccn,memdMMh2o,memdNNccn) 5 6 ! to use 'getin' 6 7 USE comcstfi_h … … 31 32 c values which are then used by the sedimentation and advection 32 33 c schemes. 34 c CO2 ice particles can nucleate on both dust and on water ice particles 35 c When CO2 ice is deposited onto a water ice particles, the particle is 36 c removed from the water tracers. 37 cWARNING: no sedimentation of the water ice origin is performed 38 c in the microphysical timestep in co2cloud.F. 33 39 34 40 c Authors of the water ice clouds microphysics … … 66 72 REAL tauscaling(ngrid) ! Convertion factor for qdust and Ndust 67 73 68 realrice(ngrid,nlay) ! Water Ice mass mean radius (m)74 REAL rice(ngrid,nlay) ! Water Ice mass mean radius (m) 69 75 ! used for nucleation of CO2 on ice-coated ccns 70 76 … … 98 104 REAL l0,l1,l2,l3,l4 99 105 REAL cste 100 REALdMice ! mass of condensed ice106 DOUBLE PRECISION dMice ! mass of condensed ice 101 107 DOUBLE PRECISION sumcheck 102 108 DOUBLE PRECISION pco2 ! Co2 vapor partial pressure (Pa) … … 104 110 DOUBLE PRECISION Mo,No 105 111 DOUBLE PRECISION Rn, Rm, dev2, n_derf, m_derf 106 112 DOUBLE PRECISION memdMMccn(ngrid,nlay) 113 DOUBLE PRECISION memdMMh2o(ngrid,nlay) 114 DOUBLE PRECISION memdNNccn(ngrid,nlay) 115 107 116 ! Radius used by the microphysical scheme (m) 108 117 DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin … … 116 125 c EXTERNAL sigco2 117 126 118 DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM 127 DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o 119 128 DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate 120 129 DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate 121 130 REAL seq 122 131 DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) 123 132 DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m) 124 133 REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3) 134 125 135 REAL rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3) 126 136 REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m) … … 128 138 c REAL res ! Resistance growth 129 139 DOUBLE PRECISION Ic_rice ! Mass transfer rate CO2 ice crystal 130 131 140 DOUBLE PRECISION ratioh2o_ccn 141 DOUBLE PRECISION vo2co2 132 142 c Parameters of the size discretization 133 143 c used by the microphysical scheme 134 DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e- 9! Minimum radius (m)135 DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e- 5! Maximum radius (m)136 DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.099e-9144 DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-10 ! Minimum radius (m) 145 DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-4 ! Maximum radius (m) 146 DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-11 137 147 ! Minimum boundary radius (m) 138 DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e- 2! Maximum boundary radius (m)148 DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-3 ! Maximum boundary radius (m) 139 149 DOUBLE PRECISION vrat_cld ! Volume ratio 140 150 DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m) … … 144 154 DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3) 145 155 146 DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff 156 DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o 147 157 REAL sigma_iceco2 ! Variance of the ice and CCN distributions 148 158 SAVE sigma_iceco2 … … 223 233 c mteta = 0.952 224 234 c mtetaco2 = 0.952 225 cwrite(*,*) 'co2_param contact parameter:', mtetaco2235 write(*,*) 'co2_param contact parameter:', mtetaco2 226 236 227 237 c Volume of a co2 molecule (m3) … … 235 245 c write(*,*) 'nuice for sedimentation:', nuiceco2_sed 236 246 c write(*,*) 'Volume of a co2 molecule:', vo1co2 247 248 write(*,*) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2 249 write(*,*) 'nuice for co2 ice sedimentation:', nuiceco2_sed 250 write(*,*) 'Volume of a co2 molecule:', vo1co2 251 237 252 238 253 … … 249 264 250 265 res_out(:,:) = 0 251 rice(:,:) = 1.e- 11266 rice(:,:) = 1.e-8 252 267 riceco2(:,:) = 1.e-11 253 268 … … 261 276 & pt(1:ngrid,1:nlay) + 262 277 & pdt(1:ngrid,1:nlay) * ptimestep 263 278 c call WRITEDIAGFI(ngrid,"Ztclouds", 279 c & "Ztclouds",'K',3,zt) 280 c call WRITEDIAGFI(ngrid,"pdtclouds", 281 c & "pdtclouds",'K',3,pdt*ptimestep) 282 264 283 zq(1:ngrid,1:nlay,1:nq) = 265 284 & pq(1:ngrid,1:nlay,1:nq) + … … 286 305 287 306 c Main loop over the GCM's grid 288 DO l=1,nlay289 DO ig=1,ngrid307 DO l=1,nlay 308 DO ig=1,ngrid 290 309 291 310 … … 297 316 ! 3. Nucleation 298 317 !============================================================= 299 318 rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*zt(ig,l) 319 & -2.87e-6*zt(ig,l)*zt(ig,l)) 320 vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l)) ! m0co2 321 rho_ice_co2=rho_ice_co2T(ig,l) 300 322 c call updaterccn(zq(ig,l,igcm_dust_mass), 301 323 c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) 302 324 c write(*,*) "l, pco2, satu= ",l,pco2,satu 303 325 IF ( satu .ge. 1d0 ) THEN ! if there is condensation 304 326 c write(*,*) 305 c write(*,*) "l, pco2, satu= ",l,pco2,satu 327 328 !write(*,*) "Zt, rho=",zt(ig,l),rho_ice_co2 306 329 c Masse_atm=mmean(ig,l)*1.e-3*pplay(ig,l)/rgp/zt(ig,l) !Kg par couche 307 330 308 331 309 call updaterccn(zq(ig,l,igcm_dust_mass), 310 & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) 332 c call updaterccn(zq(ig,l,igcm_dust_mass), 333 c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) 334 335 c call updaterccn(zq(ig,l,igcm_dust_mass), 336 c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) 311 337 c write(*,*) "Improved, l,Rdust = ",l,rdust(ig,l) 312 338 … … 315 341 & / zq(ig,l,igcm_dust_number) 316 342 rdust(ig,l)= rdust(ig,l)**(1./3.) 317 cwrite(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l)318 rdust(ig,l)=max(1.e- 9,rdust(ig,l))319 rdust(ig,l)=min( 2e-6,rdust(ig,l))320 c write(*,*) "Improved3, l,Rdust = ",l,rdust(ig,l)321 343 !write(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l) 344 rdust(ig,l)=max(1.e-10,rdust(ig,l)) 345 rdust(ig,l)=min(5.e-5,rdust(ig,l)) 346 ! write(*,*) "Improved3,Rdust = ",rdust(ig,l) 347 322 348 c Expand the dust moments into a binned distribution 323 Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) 324 No = zq(ig,l,igcm_dust_number)* tauscaling(ig) 325 c write(*,*) " dust number, mass = ",349 Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30 350 No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30 351 c write(*,*) "Improved dust number, mass = ", 326 352 c & zq(ig,l,igcm_dust_number)* tauscaling(ig), 327 353 c & zq(ig,l,igcm_dust_mass)* tauscaling(ig) … … 340 366 n_aer(i) = n_aer(i) + 0.5 * No * n_derf 341 367 m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf 342 c write(*,*) "i, r b_cldco2(i) = ",i, rb_cldco2(i),n_aer(i)343 368 c write(*,*) "i, rad_cldco2(i) = ",i, rad_cldco2(i), 369 c & n_aer(i) 344 370 enddo 345 371 … … 353 379 print*, "WARNING, No sumcheck PROBLEM" 354 380 print*, "sumcheck, No",sumcheck, No 381 print*, "rdust =",rdust(ig,l) 355 382 print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l 356 383 print*, "Dust binned distribution", n_aer … … 372 399 endif 373 400 374 c Expand the water ice moments into a binned distribution 375 c For now the radius grid's bound are same as for dust 376 c min=100 nm and max=10microns 377 c might need a change if rice (water) is large (but how large?) - listo 378 379 Mo = zq(ig,l,igcm_h2o_ice)* tauscaling(ig) + 1.e-30 401 call updaterice_micro( 402 & zq(ig,l,igcm_h2o_ice), ! ice mass 403 & zq(ig,l,igcm_ccn_mass), ! ccn mass 404 & zq(ig,l,igcm_ccn_number), ! ccn number 405 & tauscaling(ig),rice(ig,l),rhocloud(ig,l)) 406 ! rice radius of CCN + H20 crystal 407 !write(*,*) "Improved1 Rice=",rice(ig,l) 408 rice(ig,l)=max(1.e-10,rice(ig,l)) 409 rice(ig,l)=min(5.e-5,rice(ig,l)) 410 !write(*,*) "Improved2 Rice=",rice(ig,l) 411 Mo = zq(ig,l,igcm_h2o_ice) + 412 & zq(ig,l,igcm_ccn_mass)*tauscaling(ig) 413 & + 1.e-30 !Total mass of H20 crystals,CCN included 380 414 No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30 381 415 Rn = rice(ig,l) … … 389 423 n_derf = derf( (rb_cldco2(i+1)+Rn) *dev2) 390 424 m_derf = derf( (rb_cldco2(i+1)+Rm) *dev2) 391 n_aer_h2oice(i) = n_aer(i) + 0.5 * No * n_derf ! vector not really needed - temp var - listo 392 m_aer_h2oice(i) = m_aer(i) + 0.5 * Mo * m_derf ! vector not really needed - temp var 393 rad_h2oice(i) = ((m_aer_h2oice(i)/rho_ice/ 394 & n_aer_h2oice(i) + vol_cld(i))) 395 & *0.75/pi**(1./3) 396 c write(*,*) "before nuc, i,rad_h2o(i)= ",i,rad_h2oice(i) 397 c & ,m_aer(i),n_aer(i) 425 n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf ! vector not really needed - temp var - listo 426 m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf ! vector not really needed - temp var 427 rad_h2oice(i) = rad_cldco2(i) 428 c write(*,*) "before nuc, i,rad_h2o(i)= ",i,rad_cldco2(i) 429 c & ,m_aer_h2oice(i),n_aer_h2oice(i) 398 430 enddo 399 431 sumcheck = 0 432 do i = 1, nbinco2_cld 433 sumcheck = sumcheck + n_aer_h2oice(i) 434 enddo 435 sumcheck = abs(sumcheck/No - 1) 436 if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then 437 print*, "WARNING, Noh2o sumcheck PROBLEM" 438 print*, "sumcheck, No",sumcheck, No 439 print*, "rice =",rice(ig,l) 440 print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l 441 print*, "Dust binned distribution", n_aer_h2oice 442 STOP 443 endif 400 444 401 445 … … 403 447 call nucleaCO2(dble(pco2),zt(ig,l),dble(satu) 404 448 & ,n_aer,rate,n_aer_h2oice 405 & ,rad_h2oice,rateh2o )449 & ,rad_h2oice,rateh2o,vo2co2) 406 450 ! regarder rateh20, et mettre = 0 si non nul pour le moment 407 451 dN = 0. … … 410 454 dMh2o = 0. 411 455 do i = 1, nbinco2_cld 412 !n_aer(i) = n_aer(i)/( 1. + (rate(i)+rateh2o(i))*ptimestep)413 !m_aer(i) = m_aer(i)/( 1. + (rate(i)+rateh2o(i))*ptimestep)414 456 Proba=1.0-dexp(-1.*ptimestep*rate(i)) 415 416 417 c dNh2o = dNh2o + n_aer_h2oice(i) * rateh2o(i) 418 c & * ptimestep 419 c dMh2o = dMh2o + m_aer_h2oice(i) * rateh2o(i) 420 c & *ptimestep 457 Probah2o=1.0-dexp(-1.*ptimestep*rateh2o(i)) 458 459 dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o 460 dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o 421 461 422 462 dN = dN + n_aer(i) * Proba … … 425 465 enddo 426 466 427 c do i=1,nbinco2_cld428 c write(*,*) "i n_aer m_aer = ",i,n_aer(i),m_aer(i)429 c enddo430 467 ! dM masse activée (kg) et dN nb particules par kg d'air 431 468 432 c Update Dust particles433 c For CO2 ice : no subtraction from dust (neither for water ice particles)434 ! zq(ig,l,igcm_dust_mass) =435 ! & zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)436 ! zq(ig,l,igcm_dust_number) =437 ! & zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)438 469 c write(*,*) " nuclea dM = ",dM/tauscaling(ig), 439 470 c & " nuclea dN = ", dN/tauscaling(ig) … … 441 472 dNN= dN/tauscaling(ig) 442 473 dMM= dM/tauscaling(ig) 443 444 dNN=min(dNN,abs(zq0(ig,l,igcm_dust_number))) 445 dMM=min(dMM,zq0(ig,l,igcm_dust_mass)) 474 dNNh2o=dNh2o/tauscaling(ig) 475 dMMh2o=dMh2o/tauscaling(ig) 476 477 dNN=min(dNN,abs(zq(ig,l,igcm_dust_number))) 478 dMM=min(dMM,abs(zq(ig,l,igcm_dust_mass))) 479 c 480 c write(*,*) "Nuclea dNN crees=",dNN 481 dNNh2o=min(dNNh2o,abs(zq(ig,l,igcm_ccn_number))) 482 dMMh2o=min(dMMh2o,abs(zq(ig,l,igcm_h2o_ice)/tauscaling(ig) 483 & +zq(ig,l,igcm_ccn_mass))) !Total mass of H2O crystals available 484 485 c write(*,*) "Nuclea dNNh2o crees=",dNNh2o 446 486 447 487 c Update CCNs for CO2 crystals … … 449 489 zq(ig,l,igcm_ccnco2_mass) = 450 490 & zq(ig,l,igcm_ccnco2_mass) + dMM 451 452 491 zq(ig,l,igcm_ccnco2_number) = 453 492 & zq(ig,l,igcm_ccnco2_number) + dNN 454 455 493 zq(ig,l,igcm_dust_mass) = 456 494 & zq(ig,l,igcm_dust_mass) - dMM … … 458 496 & zq(ig,l,igcm_dust_number) - dNN 459 497 460 461 c + enlever les CCN a la distri de dust 462 463 c write(*,*) "new dust_mass, number =", 464 c & zq(ig,l,igcm_dust_mass)* tauscaling(ig), 465 c & zq(ig,l,igcm_dust_number)*tauscaling(ig) 466 c write(*,*) "new ccn mass, number =", 467 c & zq(ig,l,igcm_ccnco2_mass)* tauscaling(ig) 468 c & ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) 498 c Update CCN for CO2 nucleating on H2O CCN : 499 ! Warning: must keep memory of it 500 zq(ig,l,igcm_ccnco2_mass) = 501 & zq(ig,l,igcm_ccnco2_mass) + dMMh2o 502 zq(ig,l,igcm_ccnco2_number) = 503 & zq(ig,l,igcm_ccnco2_number) + dNNh2o 504 505 506 zq(ig,l,igcm_ccn_number) = 507 & zq(ig,l,igcm_ccn_number) - dNNh2o 508 509 ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice) 510 & +zq(ig,l,igcm_ccn_mass)*tauscaling(ig)) 511 512 513 memdMMh2o(ig,l)= memdMMh2o(ig,l)+zq(ig,l,igcm_h2o_ice)* 514 & dMMh2o*ratioh2o_ccn 515 memdMMccn(ig,l)= memdMMccn(ig,l)+zq(ig,l,igcm_ccn_mass)* 516 & tauscaling(ig)*dMMh2o*ratioh2o_ccn 517 memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o 518 c if (dMMh2o .gt. 0) then 519 c write(*,*) 'test h2o' 520 c write(*,*) "dMMh2o=",dMMh2o 521 c write(*,*) "2 =",zq(ig,l,igcm_ccn_mass)*tauscaling(ig)* 522 c & dMMh2o*ratioh2o_ccn+zq(ig,l,igcm_h2o_ice)* 523 c & dMMh2o*ratioh2o_ccn 524 c write(*,*) "3=",zq(ig,l,igcm_ccn_mass)*tauscaling(ig)* 525 c & dMMh2o*ratioh2o_ccn 526 c write(*,*) "4=",zq(ig,l,igcm_h2o_ice)* 527 c & dMMh2o*ratioh2o_ccn 528 c write(*,*) "tauscaling=",tauscaling(ig) 529 c endif 530 zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)* 531 & (1.-dMMh2o*ratioh2o_ccn) 532 zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass)* 533 & tauscaling(ig)*(1.-dMMh2o*ratioh2o_ccn) 534 469 535 470 536 ENDIF ! of is satu >1 … … 478 544 c IF ( zq(ig,l,igcm_ccnco2_number)*tauscaling(ig).ge. 1.0) THEN 479 545 480 IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. threshJA) 481 & THEN ! we trigger crystal growth 482 483 call updaterccn(zq(ig,l,igcm_dust_mass), 484 & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) 546 IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. 1.) 547 & THEN 548 ! we trigger crystal growth 549 c 550 551 c Niceco2 = zq(ig,l,igcm_co2_ice) 552 c Qccnco2 = zq(ig,l,igcm_ccnco2_mass) 553 c Nccnco2 = zq(ig,l,igcm_ccnco2_number) 554 c call updaterice_microco2(Niceco2,Qccnco2,Nccnco2, 555 c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 556 c write(*,*) "updater rice=",riceco2(ig,l) 557 485 558 rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass) 486 559 & *0.75/pi/rho_dust 487 560 & / zq(ig,l,igcm_ccnco2_number) 488 rdust(ig,l)= rdust(ig,l)**(1./3.) 489 561 rdust(ig,l)= rdust(ig,l)**(1./3.) 490 562 rdust(ig,l)=max(1.e-10,rdust(ig,l)) 491 rdust(ig,l)=min(2e-6,rdust(ig,l)) 492 ! rdust(ig,l)=1.e-7 493 494 IF (zq(ig,l,igcm_ccnco2_mass) .lt. 0. .or. 495 & zq(ig,l,igcm_ccnco2_number) .lt. 0. .or. 496 & zq(ig,l,igcm_dust_mass) .lt. 0. .or. 497 & zq(ig,l,igcm_dust_number) .lt. 0. ) THEN 498 499 write(*,*) "CO2 clouds before growth CCN N,M = " 500 & ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) 501 & ,zq(ig,l,igcm_ccnco2_mass)*tauscaling(ig) 502 503 write(*,*) "before growth dust number mass = ", 504 & zq(ig,l,igcm_dust_number)*tauscaling(ig), 505 & zq(ig,l,igcm_dust_mass)*tauscaling(ig) 506 STOP 507 END IF 508 509 c write(*,*) "reff dN = ",reff,dN 510 c reff=reff/dble(dN) 511 c if (zq(ig,l,igcm_co2_ice) .le. 1.e-20) then 512 c riceco2(ig,l)=reff 513 c endif 514 515 c write(*,*) "Rdust in improved = ",rdust(ig,l) 516 563 ! rdust(ig,l)=min(5.e-6,rdust(ig,l)) 564 517 565 riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ 518 & (4.0*rho_ice_co2* pi*zq(ig,l,igcm_ccnco2_number)519 & * tauscaling(ig)) +rdust(ig,l)*rdust(ig,l)566 & (4.0*rho_ice_co2*zq(ig,l,igcm_ccnco2_number) 567 & *pi*tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) 520 568 & *rdust(ig,l) )**(1.0/3.0) 521 569 570 c riceco2(ig,l)=max(1.e-10,riceco2(ig,l)) 571 c riceco2(ig,l)=min(1.e-5,riceco2(ig,l)) 522 572 ! WATCH OUT: CO2 nuclei is supposed to be dust 523 573 ! only when deriving rhocloud (otherwise would need to keep info on water embedded in co2) - listo … … 526 576 527 577 !! Niceco2,Qccnco2,Nccnco2 528 Niceco2 = zq(ig,l,igcm_co2_ice) 529 Qccnco2 = zq(ig,l,igcm_ccnco2_mass) 530 Nccnco2 = zq(ig,l,igcm_ccnco2_number) 531 call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, 532 & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 533 c write(*,*) "Riceco2 update before growth = ",riceco2(ig,l) 534 535 No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig) 536 & + 1.e-30 578 c Niceco2 = zq(ig,l,igcm_co2_ice) 579 c Qccnco2 = zq(ig,l,igcm_ccnco2_mass) 580 c Nccnco2 = zq(ig,l,igcm_ccnco2_number) 581 c call updaterice_microco2(Niceco2,Qccnco2,Nccnco2, 582 c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 583 !write(*,*) "Riceco2 update before growth = ",riceco2(ig,l) 584 c write(*,*) "rdust before growth = ",rdust(ig,l) 585 c write(*,*) "co2 before growth=",zq(ig,l,igcm_co2) 586 c write(*,*) "pplay before growth=",pplay(ig,l) 587 c write(*,*) "zt before growth =",zt(ig,l) 588 ! write(*,*) "Satu before growth=",satu 589 c riceco2(ig,l)=max(riceco2(ig,l),rdust(ig,l)) 590 No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30 537 591 ! No nb de particules de poussieres mis à l'échelle pour donner une opacité optique 538 592 … … 549 603 drsurdt=-1.0/(4.0*pi*riceco2(ig,l)* 550 604 & riceco2(ig,l)*rho_ice_co2)*Ic_rice 551 dMice = No * Ic_rice *ptimestep ! Kg par kg d'air, <0 si croissance !605 dMice = No * Ic_rice*ptimestep ! Kg par kg d'air, <0 si croissance ! 552 606 c write(*,*) "dMicev0 in improved = " , dMice 553 607 554 if (dMice .gt. 0) dMice = min(dMice,zq0(ig,l,igcm_co2_ice)) 555 if (dMice .lt. 0) dMice = max(dMice,-1.*zq0(ig,l,igcm_co2)) 556 557 558 559 560 riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep 608 if (dMice .ge. 0d0) then 609 dMice = min(dMice,abs(zq(ig,l,igcm_co2_ice))) 610 else 611 dMice =-1.* min(abs(dMice),abs(zq(ig,l,igcm_co2))) 612 endif 613 riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep 561 614 c write(*,*) "riceco2+dr/dt = ", riceco2(ig,l) 562 563 c write(*,*) "dMice in improved = " , dMice 615 c write(*,*) "dMice in improved = " , dMice 564 616 565 566 zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice) 567 & -dMice 568 zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)+dMice 569 c write(*,*) "zq co2 ice = ", zq(ig,l,igcm_co2_ice) 570 countcells = countcells + 1 617 zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice) 618 & -dMice 619 zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)+dMice 620 c write(*,*) "Improved zq co2 ice = ", zq(ig,l,igcm_co2_ice) 621 ! countcells = countcells + 1 571 622 572 riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ 573 & (4.0*rho_ice_co2*pi*zq(ig,l,igcm_ccnco2_number) 574 & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) 575 & *rdust(ig,l) )**(1.0/3.0) 576 c write(*,*) "new riceco2 = ",riceco2(ig,l) 577 578 !! Niceco2,Qccnco2,Nccnco2 579 Niceco2 = zq(ig,l,igcm_co2_ice) 580 Qccnco2 = zq(ig,l,igcm_ccnco2_mass) 581 Nccnco2 = zq(ig,l,igcm_ccnco2_number) 582 call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, 583 & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 584 c write(*,*) "new riceco2 updaterad= ",riceco2(ig,l) 623 c riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ 624 c & (4.0*rho_ice_co2*pi*zq(ig,l,igcm_ccnco2_number) 625 c & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) 626 c & *rdust(ig,l) )**(1.0/3.0) 627 c write(*,*) "Improved new riceco2 = ",riceco2(ig,l) 628 629 c write(*,*) "new riceco2 improvedupdaterad= ",riceco2(ig,l) 585 630 586 631 ! latent heat release … … 598 643 pdtcloudco2(ig,l)= dMice*lw/cpp/ptimestep ! kg par couche * J par kg /J par K / s = K par seconde 599 644 600 c write(*,*) "pdtcloudco2 = ",pdtcloudco2(ig,l)645 c write(*,*) "pdtcloudco2 after growth = ",pdtcloudco2(ig,l) 601 646 602 c Peut etre -1*dMice?603 604 647 605 648 !deltaT par condens/subli. qui remplace le dT du CO2 du newcondens pré-constantino … … 618 661 619 662 c interaction ho-co2 ici, dans la mise a jour des traceurs WARNING reflechir 663 !! Niceco2,Qccnco2,Nccnco2 664 665 620 666 621 622 ENDIF !of if Nccn>1 623 c if (riceco2(ig,l) .lt. 1.e-9) then 624 625 if (zq(ig,l,igcm_co2_ice).le.1.e-20 .or. 626 & riceco2(ig,l) .lt. 1.e-9 .or. riceco2(ig,l) 627 & .ge. 4.999e-4) then 628 ! Reverser les ccn libérés dans les h2o ou dust? 629 630 c c ICI: Co2 ice devient vapeur 631 zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) 632 & + zq(ig,l,igcm_co2_ice) 633 634 zq(ig,l,igcm_dust_mass) = 635 & zq(ig,l,igcm_dust_mass) 636 & + zq(ig,l,igcm_ccnco2_mass) 637 zq(ig,l,igcm_dust_number) 638 & = zq(ig,l,igcm_dust_number) 639 & + zq(ig,l,igcm_ccnco2_number) 640 c c CCNs 641 zq(ig,l,igcm_ccnco2_mass) = 0. 642 zq(ig,l,igcm_ccnco2_number) =0. 643 zq(ig,l,igcm_co2_ice) = 0. 644 riceco2(ig,l)=0. 645 endif 646 647 c write(*,*) "zq co2 end imp= ", zq(i g,l,igcm_co2_ice),satu 648 667 rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass) 668 & *0.75/pi/rho_dust 669 & / zq(ig,l,igcm_ccnco2_number) 670 rdust(ig,l)= rdust(ig,l)**(1./3.) 671 rdust(ig,l)=max(1.e-10,rdust(ig,l)) 672 !rdust(ig,l)=min(5.e-6,rdust(ig,l)) 673 674 riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ 675 & (4.0*rho_ice_co2*zq(ig,l,igcm_ccnco2_number) 676 & *pi*tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) 677 & *rdust(ig,l) )**(1.0/3.0) 678 !Niceco2 = zq(ig,l,igcm_co2_ice) 679 !Qccnco2 = zq(ig,l,igcm_ccnco2_mass) 680 !Nccnco2 = zq(ig,l,igcm_ccnco2_number) 681 c 682 c call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, 683 c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 649 684 650 651 ENDDO ! of ig loop 652 ENDDO ! of nlayer loop 685 if ((zq(ig,l,igcm_co2_ice).le. 1.e-23) 686 & .or.(riceco2(ig,l) .le. rdust(ig,l))) then 687 688 c write(*,*) "Riceco2 improved before reset=",riceco2(ig,l) 689 c write(*,*) "Niceco2 improved before reset=", 690 c & zq(ig,l,igcm_co2_ice) 691 c write(*,*) "Rdust improved before reset=",rdust(ig,l) 692 693 if (memdMMccn(ig,l) .gt. 0) then 694 zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass) 695 & +memdMMccn(ig,l) 696 endif 697 if (memdMMh2o(ig,l) .gt. 0) then 698 zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice) 699 & +memdMMh2o(ig,l) 700 endif 701 702 if (memdNNccn(ig,l) .gt. 0) then 703 zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number) 704 & +memdNNccn(ig,l) 705 endif 706 707 if (zq(ig,l,igcm_ccnco2_mass) .gt. 1.e-30) then 708 zq(ig,l,igcm_dust_mass) = 709 & zq(ig,l,igcm_dust_mass) 710 & + zq(ig,l,igcm_ccnco2_mass)- 711 & (memdMMh2o(ig,l)+memdMMccn(ig,l)) 712 endif 713 if (zq(ig,l,igcm_ccnco2_number) .gt. 1.e-30) then 714 zq(ig,l,igcm_dust_number) = 715 & zq(ig,l,igcm_dust_number) 716 & + zq(ig,l,igcm_ccnco2_number)-memdNNccn(ig,l) 717 endif 718 719 if (zq(ig,l,igcm_co2_ice) .gt. 1.e-30) then 720 zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) 721 & + zq(ig,l,igcm_co2_ice) 722 endif 723 724 zq(ig,l,igcm_ccnco2_mass)=0. 725 zq(ig,l,igcm_co2_ice)=0. 726 zq(ig,l,igcm_ccnco2_number)=0. 727 memdNNccn(ig,l)=0. 728 memdMMh2o(ig,l)=0. 729 memdMMccn(ig,l)=0. 730 riceco2(ig,l)=0. 731 pdtcloudco2(ig,l)=0. 732 endif 733 734 ENDIF ! of if NCCN > 1 735 ENDDO ! of ig loop 736 ENDDO ! of nlayer loop 653 737 654 738 … … 661 745 & (zq(1:ngrid,1:nlay,igcm_co2_ice) - 662 746 & zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep 663 c do l=1,nlay 664 c write(*,*) "end imp",pdqcloudco2(1,l,igcm_co2), 665 c & pdqcloudco2(1,l,igcm_co2_ice) 666 c enddo 747 748 pdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) = 749 & (zq(1:ngrid,1:nlay,igcm_h2o_ice) - 750 & zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep 751 752 pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) = 753 & (zq(1:ngrid,1:nlay,igcm_ccn_mass) - 754 & zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep 755 756 pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) = 757 & (zq(1:ngrid,1:nlay,igcm_ccn_number) - 758 & zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep 759 667 760 pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) = 668 761 & (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) - … … 672 765 & (zq(1:ngrid,1:nlay,igcm_ccnco2_number) - 673 766 & zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep 674 767 675 768 676 if (scavenging) then769 c if (scavenging) then 677 770 678 771 pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) = 679 772 & (zq(1:ngrid,1:nlay,igcm_dust_mass) - 680 773 & zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep 774 681 775 pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) = 682 776 & (zq(1:ngrid,1:nlay,igcm_dust_number) - 683 777 & zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep 684 endif 685 686 c call WRITEdiagfi(ngrid,"Improvedb","co2 traceur","kg/kg",1, 687 c & zq0(1,:,igcm_co2_ice) ) 688 689 c call WRITEdiagfi(ngrid,"Improveda","co2 traceur","kg/kg",1, 690 c & zq(1,:,igcm_co2_ice) ) 691 692 693 694 695 c call WRITEdiagfi(ngrid,"satuco2","satu co2 in improved","kg/kg",1, 696 c & satu) 697 698 699 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 700 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 701 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 702 IF (test_flag) then 703 704 ! error2d(:) = 0. 705 DO l=1,nlay 706 DO ig=1,ngrid 707 ! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig)) 708 satubf(ig,l) = zq0(ig,l,igcm_co2)/zqsat(ig,l) ! att. if zqsat=mmr or psat 709 satuaf(ig,l) = zq(ig,l,igcm_co2)/zqsat(ig,l) 710 ENDDO 711 ENDDO 712 713 write(*,*) 'count is ',countcells, ' i.e. ', 714 & countcells*100/(nlay*ngrid), '% for microphys computation' 715 716 ! IF (ngrid.ne.1) THEN ! 3D 717 ! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3, 718 ! & satu_out) 719 ! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3, 720 ! & dM_out) 721 ! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3, 722 ! & dN_out) 723 ! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2, 724 ! & error2d) 725 ! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3, 726 ! & zqsat) 727 ! ENDIF 728 729 ! IF (ngrid.eq.1) THEN ! 1D 730 ! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1, 731 ! & error_out) 732 ! call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1, 733 ! & res_out) 734 call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1, 735 & satubf) 736 call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1, 737 & satuaf) 738 call WRITEdiagfi(ngrid,"vapbf","CO2vap before","kg/kg",1, 739 & zq0(1,1,igcm_co2)) 740 call WRITEdiagfi(ngrid,"vapaf","CO2vap after","kg/kg",1, 741 & zq(1,1,igcm_co2)) 742 call WRITEdiagfi(ngrid,"icebf","CO2ice before","kg/kg",1, 743 & zq0(1,1,igcm_co2_ice)) 744 call WRITEdiagfi(ngrid,"iceaf","CO2ice after","kg/kg",1, 745 & zq(1,1,igcm_co2_ice)) 746 call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1, 747 & zq0(1,1,igcm_ccnco2_number)) 748 call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1, 749 & zq(1,1,igcm_ccnco2_number)) 750 c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1, 751 c & gr_out) 752 c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1, 753 c & rate_out) 754 c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1, 755 c & dM_out) 756 c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1, 757 c & dN_out) 758 c call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1, 759 c & zqsat) 760 ! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1, 761 ! & satu_out) 762 ! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1, 763 ! & rdust) 764 ! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1, 765 ! & rsedcloud) 766 ! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1, 767 ! & rhocloud) 768 ! ENDIF 769 770 ENDIF ! endif test_flag 771 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 772 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 773 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 774 778 c endif 779 775 780 return 776 781 end … … 778 783 779 784 780 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc781 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc782 c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0783 c It is an analytical solution to the ice radius growth equation,784 c with the approximation of a constant 'reduced' cunningham correction factor785 c (lambda in growthrate.F) taken at radius req instead of rice786 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc787 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc788 789 c subroutine phi(rice,req,coeff1,coeff2,time)790 c791 c implicit none792 c793 c ! inputs794 c real rice ! ice radius795 c real req ! ice radius at equilibirum796 c real coeff1 ! coeff for the log797 c real coeff2 ! coeff for the arctan798 c799 c ! output800 c real time801 c802 c !local803 c real var804 c805 c ! 1.73205 is sqrt(3)806 c807 c var = max(808 c & abs(rice-req) / sqrt(rice*rice + rice*req + req*req),1e-30)809 c810 c time =811 c & coeff1 *812 c & log( var )813 c & + coeff2 * 1.73205 *814 c & atan( (2*rice+req) / (1.73205*req) )815 c816 c return817 c end818 819 820 -
trunk/LMDZ.MARS/libf/phymars/initracer.F
r1660 r1685 375 375 rho_dust=2500. ! Mars dust density (kg.m-3) 376 376 rho_ice=920. ! Water ice density (kg.m-3) 377 rho_ice_co2=1500. !dry ice density (kg.m-3), varies with T from 0.98 to 1.5 see Satorre et al., PSS 2008 377 rho_ice_co2=1650. 378 !Mangan et al., Icarus 2017 :CO2 density = 1.72391-2.53×10−4T – 2.87×10−6T^2 378 379 nuice_ref=0.1 ! Effective variance nueff of the 379 380 ! water-ice size distribution … … 527 528 alpha_lift(igcm_co2) =0. 528 529 alpha_devil(igcm_co2)=0. 529 radius(igcm_co2_ice)=1.e- 6530 radius(igcm_co2_ice)=1.e-8 530 531 rho_q(igcm_co2_ice)=rho_ice_co2 531 532 alpha_lift(igcm_co2_ice) =0. -
trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F
r1620 r1685 52 52 53 53 54 Tcm = dble(T) ! initialize pourquoi 0 et pas t(i)54 Tcm = 0!dble(T) ! initialize pourquoi 0 et pas t(i) 55 55 56 56 T_inf = 0d0 57 T_sup = 200d057 T_sup = 800d0 58 58 59 59 T_dT = 0.1 ! precision - mettre petit et limiter nb iteration? -
trunk/LMDZ.MARS/libf/phymars/microphys.h
r1617 r1685 22 22 DOUBLE PRECISION, PARAMETER :: mn2 = 28.01d-3 23 23 ! Effective CO2 gas molecular radius (m) 24 DOUBLE PRECISION, PARAMETER :: molco2 = 2.2d-10 24 ! bachnar 2016 value :1.97d-10 ! old value = 2.2d-10 25 DOUBLE PRECISION, PARAMETER :: molco2 = 1.97d-10 25 26 ! Effective H2O gas molecular radius (m) 26 27 DOUBLE PRECISION, PARAMETER :: molh2o = 1.2d-10 … … 54 55 !CO2 part 55 56 ! number of bins for nucleation 56 INTEGER, PARAMETER :: nbinco2_cld= 4057 INTEGER, PARAMETER :: nbinco2_cld=100 57 58 ! Surface tension of ice/vapor (J.m-2) 58 59 DOUBLE PRECISION, PARAMETER :: sigco2 = 0.08 … … 60 61 ! water on a dust-like substrate 61 62 ! (J/molecule) 62 DOUBLE PRECISION, PARAMETER :: desorpco2 = 3.25e-20 63 DOUBLE PRECISION, PARAMETER ::desorpco2=3.07e-20 64 ! bachnar 2016 value 3.07d-20 65 !old value 3.20e-20 63 66 ! Jump frequency of a co2 molecule (s-1) 64 67 DOUBLE PRECISION, PARAMETER :: nusco2 = 2.9e+12 … … 71 74 ! Contact parameter ( m=cos(theta) ) 72 75 ! (initialized in improvedCO2clouds.F) 73 REAL, parameter :: mtetaco2 = 0.952 76 ! bachnar 2016 value :0.78 77 !old value 0.952 78 REAL, parameter :: mtetaco2 = 0.78 74 79 ! Volume of a co2 molecule (m3) 75 80 DOUBLE PRECISION :: vo1co2 76 81 ! Radius used by the microphysical scheme (m) 77 82 DOUBLE PRECISION :: rad_cldco2(nbinco2_cld) 78 REAL, parameter :: threshJA = 1 .083 REAL, parameter :: threshJA = 1 79 84 ! COMMON/microphys/vo1co2,rad_cldco2 80 85 -
trunk/LMDZ.MARS/libf/phymars/nucleaCO2.F
r1617 r1685 2 2 * * 3 3 subroutine nucleaCO2(pco2,temp,sat,n_ccn,nucrate, 4 & n_ccn_h2oice,rad_h2oice,nucrate_h2oice) 4 & n_ccn_h2oice,rad_h2oice,nucrate_h2oice, 5 & vo2co2) 5 6 USE comcstfi_h 6 7 … … 24 25 25 26 c Inputs 26 DOUBLE PRECISION pco2,sat 27 DOUBLE PRECISION pco2,sat,vo2co2 27 28 DOUBLE PRECISION n_ccn(nbinco2_cld), n_ccn_h2oice(nbinco2_cld) 28 29 REAL temp … … 49 50 c double precision xratio 50 51 51 double precision mtetalocal ! local mteta in double precision52 double precision mtetalocal,mtetalocalh ! local mteta in double precision 52 53 53 54 double precision fshapeco2simple,zefshapeco2 … … 63 64 64 65 mtetalocal = dble(mtetaco2) !! use mtetalocal for better performance 66 mtetalocalh=dble(mteta) 65 67 66 68 cccccccccccccccccccccccccccccccccccccccccccccccccc … … 99 101 100 102 nco2 = pco2 / kbz / temp 101 rstar = 2. * sigco2 * vo 1co2 / (kbz*temp*dlog(sat))102 gstar = 4. * pi * (rstar * rstar * rstar) / (3.*vo 1co2)103 rstar = 2. * sigco2 * vo2co2 / (kbz*temp*dlog(sat)) 104 gstar = 4. * pi * (rstar * rstar * rstar) / (3.*vo2co2) 103 105 104 106 fshapeco2simple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal) … … 152 154 153 155 if (rad_h2oice(i).gt.3000.*rstar) then 154 zefshapeco2 = fshapeco2simple 155 else 156 zefshapeco2 = fshapeco2(mtetalocal,rad_h2oice(i)/rstar) ! same m for dust/h2o ice 156 zefshapeco2 = (2.+mtetalocalh)*(1.-mtetalocalh)* 157 & (1.-mtetalocalh) / 4. 158 else 159 zefshapeco2 = fshapeco2(mtetalocalh,rad_h2oice(i)/rstar) ! same m for dust/h2o ice 157 160 endif 158 161 -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r1684 r1685 1164 1164 & pq,pdq,zdqcloudco2,zdtcloudco2, 1165 1165 & nq,tau,tauscaling,rdust,rice,riceco2,nuice, 1166 & rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2) 1166 & rsedcloudco2,rhocloudco2, 1167 & rsedcloud,rhocloud,zzlev,zdqssed_co2) 1167 1168 1168 1169
Note: See TracChangeset
for help on using the changeset viewer.
