Changeset 2953 for trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F
- Timestamp:
- Apr 28, 2023, 2:31:03 PM (19 months ago)
- File:
-
- 1 edited
-
trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F
r2942 r2953 8 8 CONTAINS 9 9 10 SUBROUTINE co2condens(ngrid,nlayer,nq, ptimestep,10 SUBROUTINE co2condens(ngrid,nlayer,nq,nslope,ptimestep, 11 11 $ pcapcal,pplay,pplev,ptsrf,pt, 12 12 $ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, … … 34 34 USE vertical_layers_mod, ONLY: ap, bp 35 35 #endif 36 use comslope_mod, ONLY: subslope_dist,def_slope_mean 36 37 IMPLICIT NONE 37 38 c======================================================================= … … 59 60 INTEGER,INTENT(IN) :: nlayer ! number of vertical layers 60 61 INTEGER,INTENT(IN) :: nq ! number of tracers 62 INTEGER,INTENT(IN) :: nslope ! number of subslope 61 63 62 64 REAL,INTENT(IN) :: ptimestep ! physics timestep (s) 63 REAL,INTENT(IN) :: pcapcal(ngrid )65 REAL,INTENT(IN) :: pcapcal(ngrid,nslope) 64 66 REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa) 65 67 REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) 66 REAL,INTENT(IN) :: ptsrf(ngrid ) ! surface temperature (K)68 REAL,INTENT(IN) :: ptsrf(ngrid,nslope) ! surface temperature (K) 67 69 REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K) 68 70 REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2) … … 73 75 REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2) 74 76 ! from previous physical processes 75 REAL,INTENT(IN) :: pdtsrf(ngrid ) ! tendency on surface temperature from77 REAL,INTENT(IN) :: pdtsrf(ngrid,nslope) ! tendency on surface temperature from 76 78 ! previous physical processes (K/s) 77 79 REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s) … … 84 86 REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1) 85 87 86 REAL,INTENT(INOUT) :: piceco2(ngrid ) ! CO2 ice on the surface (kg.m-2)87 REAL,INTENT(INOUT) :: psolaralb(ngrid,2 ) ! albedo of the surface88 REAL,INTENT(INOUT) :: pemisurf(ngrid ) ! emissivity of the surface88 REAL,INTENT(INOUT) :: piceco2(ngrid,nslope) ! CO2 ice on the surface (kg.m-2) 89 REAL,INTENT(INOUT) :: psolaralb(ngrid,2,nslope) ! albedo of the surface 90 REAL,INTENT(INOUT) :: pemisurf(ngrid,nslope) ! emissivity of the surface 89 91 REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius 90 92 91 93 ! tendencies due to CO2 condensation/sublimation: 92 94 REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s) 93 REAL,INTENT(OUT) :: pdtsrfc(ngrid ) ! tendency on surface temperature (K/s)95 REAL,INTENT(OUT) :: pdtsrfc(ngrid,nslope) ! tendency on surface temperature (K/s) 94 96 REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s) 95 97 REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2) … … 112 114 REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm) 113 115 REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface) 114 REAL zdiceco2(ngrid) 116 REAL zdiceco2(ngrid,nslope) 117 REAL zdiceco2_mesh_avg(ngrid) 115 118 REAL zcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s) 116 REAL zcondices(ngrid) ! condensation rate on the ground (kg/m2/s) 119 REAL zcondices(ngrid,nslope) ! condensation rate on the ground (kg/m2/s) 120 REAL zcondices_mesh_avg(ngrid) ! condensation rate on the ground (kg/m2/s) 117 121 REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s) 118 122 REAL condens_layer(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s) … … 122 126 REAL zu(nlayer),zv(nlayer) 123 127 REAL zqc(nlayer,nq),zq1(nlayer) 124 REAL ztsrf(ngrid )128 REAL ztsrf(ngrid,nslope) 125 129 REAL ztc(nlayer), ztm(nlayer+1) 126 130 REAL zum(nlayer+1) , zvm(nlayer+1) … … 129 133 REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix 130 134 REAL availco2 131 LOGICAL condsub(ngrid )132 133 real :: emisref(ngrid )135 LOGICAL condsub(ngrid,nslope) 136 137 real :: emisref(ngrid,nslope) 134 138 135 139 REAL zdq_scav(ngrid,nlayer,nq) ! tendency due to scavenging by co2 … … 170 174 REAL masseq(nlayer),wq(nlayer+1) 171 175 INTEGER ifils,iq2 176 177 c Subslope: 178 179 REAL :: emisref_tmp(ngrid) 180 REAL :: alb_tmp(ngrid,2) ! local 181 REAL :: zcondices_tmp(ngrid) ! local 182 REAL :: piceco2_tmp(ngrid) ! local 183 REAL :: pemisurf_tmp(ngrid)! local 184 LOGICAL :: condsub_tmp(ngrid) !local 185 REAL :: zfallice_tmp(ngrid,nlayer+1) 186 REAL :: condens_layer_tmp(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s) 187 INTEGER :: islope 172 188 c---------------------------------------------------------------------- 173 189 … … 231 247 pdqc(1:ngrid,1:nlayer,1:nq) = 0 232 248 233 zcondices(1:ngrid) = 0. 234 pdtsrfc(1:ngrid) = 0. 249 zcondices(1:ngrid,1:nslope) = 0. 250 zcondices_mesh_avg(1:ngrid)=0. 251 pdtsrfc(1:ngrid,1:nslope) = 0. 235 252 pdpsrf(1:ngrid) = 0. 236 condsub(1:ngrid) = .false. 237 zdiceco2(1:ngrid) = 0. 253 condsub(1:ngrid,1:nslope) = .false. 254 zdiceco2(1:ngrid,1:nslope) = 0. 255 zdiceco2_mesh_avg(1:ngrid)=0. 238 256 239 257 zfallheat=0 … … 301 319 pdtc(ig,l)=0. 302 320 IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN 303 condsub(ig )=.true.321 condsub(ig,:)=.true. 304 322 IF (zfallice(ig,l+1).gt.0) then 305 323 zfallheat=zfallice(ig,l+1)* … … 347 365 condens_column(ig) = sum(condens_layer(ig,:)) 348 366 zfallice(ig,1) = zdqssed_co2(ig) 349 piceco2(ig) = piceco2(ig) + zdqssed_co2(ig)*ptimestep 367 DO islope = 1,nslope 368 piceco2(ig,islope) = piceco2(ig,islope) + 369 & zdqssed_co2(ig)*ptimestep * 370 & cos(pi*def_slope_mean(islope)/180.) 371 ENDDO 350 372 ENDDO 351 373 call write_output("co2condens_zfallice", … … 371 393 ztcondsol(ig)= 372 394 & 1./(bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1))) 373 ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep 395 DO islope = 1,nslope 396 ztsrf(ig,islope) = ptsrf(ig,islope) + 397 & pdtsrf(ig,islope)*ptimestep 398 ENDDO 374 399 ENDDO 375 400 … … 388 413 icap=1 389 414 ENDIF 415 416 DO islope = 1,nslope 417 c Need first to put piceco2_slope(ig,islope) in kg/m^2 flat 418 419 piceco2(ig,islope) = piceco2(ig,islope) 420 & /cos(pi*def_slope_mean(islope)/180.) 421 390 422 c 391 423 c Loop on where we have condensation/ sublimation 392 IF ((ztsrf(ig ) .LT. ztcondsol(ig)) .OR. ! ground cond424 IF ((ztsrf(ig,islope) .LT. ztcondsol(ig)) .OR. ! ground cond 393 425 $ (zfallice(ig,1).NE.0.) .OR. ! falling snow 394 $ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim. 395 $ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN 396 condsub(ig) = .true. 426 $ ((ztsrf(ig,islope) .GT. ztcondsol(ig)) .AND. ! ground sublim. 427 $ ((piceco2(ig,islope)+zfallice(ig,1)*ptimestep) 428 $ .NE. 0.))) THEN 429 condsub(ig,islope) = .true. 397 430 398 431 IF (zfallice(ig,1).gt.0 .and. .not. co2clouds) then … … 406 439 c condensation or partial sublimation of CO2 ice 407 440 c """"""""""""""""""""""""""""""""""""""""""""""" 408 zcondices(ig) = pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) 409 & /(latcond*ptimestep) - zfallheat 410 pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep 441 if(ztsrf(ig,islope).LT. ztcondsol(ig)) then 442 c condensation 443 zcondices(ig,islope)=pcapcal(ig,islope) 444 & *(ztcondsol(ig)-ztsrf(ig,islope)) 445 c & /((latcond+cpp*(zt(ig,1)-ztcondsol(ig)))*ptimestep) 446 & /(latcond*ptimestep) 447 & - zfallheat 448 else 449 c sublimation 450 zcondices(ig,islope)=pcapcal(ig,islope) 451 & *(ztcondsol(ig)-ztsrf(ig,islope)) 452 & /(latcond*ptimestep) 453 & - zfallheat 454 endif 455 pdtsrfc(ig,islope) = (ztcondsol(ig) - ztsrf(ig,islope)) 456 & /ptimestep 411 457 #ifdef MESOSCALE 412 458 print*, "not enough CO2 tracer in 1st layer to condense" … … 421 467 availco2 = pq(ig,1,ico2)*((ap(1)-ap(2))+ 422 468 & (bp(1)-bp(2))*(pplev(ig,1)/g - 423 & (zcondices(ig) + zfallice(ig,1))*ptimestep)) 424 if ((zcondices(ig) + condens_layer(ig,1))*ptimestep 469 & (zcondices(ig,islope) + zfallice(ig,1)) 470 & *ptimestep)) 471 if ((zcondices(ig,islope) + condens_layer(ig,1)) 472 & *ptimestep 425 473 & .gt.availco2) then 426 zcondices(ig ) = availco2/ptimestep -474 zcondices(ig,islope) = availco2/ptimestep - 427 475 & condens_layer(ig,1) 428 pdtsrfc(ig ) = (latcond/pcapcal(ig))*429 & (zcondices(ig )+zfallheat)476 pdtsrfc(ig,islope) = (latcond/pcapcal(ig,islope))* 477 & (zcondices(ig,islope)+zfallheat) 430 478 end if 431 479 #else … … 437 485 #endif 438 486 439 c If the entire CO2 ice layer sublimes 487 c If the entire CO2 ice layer sublimes on the slope 440 488 c """""""""""""""""""""""""""""""""""""""""""""""""""" 441 489 c (including what has just condensed in the atmosphere) 442 490 if (co2clouds) then 443 IF((piceco2(ig )/ptimestep).LE.444 & -zcondices(ig ))THEN445 zcondices(ig ) = -piceco2(ig)/ptimestep446 pdtsrfc(ig )=(latcond/pcapcal(ig)) *447 & (zcondices(ig )+zfallheat)491 IF((piceco2(ig,islope)/ptimestep).LE. 492 & -zcondices(ig,islope))THEN 493 zcondices(ig,islope) = -piceco2(ig,islope)/ptimestep 494 pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) * 495 & (zcondices(ig,islope)+zfallheat) 448 496 END IF 449 497 else 450 IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE. 451 & -zcondices(ig))THEN 452 zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) 453 pdtsrfc(ig)=(latcond/pcapcal(ig)) * 454 & (zcondices(ig)+zfallheat) 498 IF((piceco2(ig,islope)/ptimestep+zfallice(ig,1)).LE. 499 & -zcondices(ig,islope))THEN 500 zcondices(ig,islope) = -piceco2(ig,islope) 501 & /ptimestep - zfallice(ig,1) 502 pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) * 503 & (zcondices(ig,islope)+zfallheat) 455 504 END IF 456 505 end if 457 506 458 c Changing CO2 ice amount and pressure :507 c Changing CO2 ice amount and pressure per slope: 459 508 c """""""""""""""""""""""""""""""""""" 460 zdiceco2(ig ) = zcondices(ig) +zfallice(ig,1)509 zdiceco2(ig,islope) = zcondices(ig,islope)+zfallice(ig,1) 461 510 & + condens_column(ig) 462 511 if (co2clouds) then 463 512 ! add here only direct condensation/sublimation 464 piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep 513 piceco2(ig,islope) = piceco2(ig,islope) + 514 & zcondices(ig,islope)*ptimestep 465 515 else 466 516 ! add here also CO2 ice in the atmosphere 467 piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep 517 piceco2(ig,islope) = piceco2(ig,islope) + 518 & zdiceco2(ig,islope)*ptimestep 468 519 end if 469 pdpsrf(ig) = -zdiceco2(ig)*g 520 521 zcondices_mesh_avg(ig) = zcondices_mesh_avg(ig) + 522 & zcondices(ig,islope)* subslope_dist(ig,islope) 523 524 zdiceco2_mesh_avg(ig) = zdiceco2_mesh_avg(ig) + 525 & zdiceco2(ig,islope)* subslope_dist(ig,islope) 526 527 END IF ! if there is condensation/sublimation 528 529 piceco2(ig,islope) = piceco2(ig,islope) 530 & * cos(pi*def_slope_mean(islope)/180.) 531 532 ENDDO !islope 533 534 pdpsrf(ig) = -zdiceco2_mesh_avg(ig)*g 470 535 471 536 IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN … … 480 545 & 'condensing more than total mass', 1) 481 546 ENDIF 482 END IF ! if there is condensation/sublimation 547 483 548 ENDDO ! of DO ig=1,ngrid 549 484 550 485 551 c ******************************************************************** … … 489 555 ! Check that amont of CO2 ice is not problematic 490 556 DO ig=1,ngrid 491 if(.not.piceco2(ig).ge.0.) THEN 492 if(piceco2(ig).le.-5.e-8) print*, 493 $ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig) 494 piceco2(ig)=0. 557 DO islope = 1,nslope 558 if(.not.piceco2(ig,islope).ge.0.) THEN 559 if(piceco2(ig,islope).le.-5.e-8) print*, 560 $ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig,islope) 561 piceco2(ig,islope)=0. 495 562 endif 563 ENDDO 496 564 ENDDO 497 565 498 566 ! Set albedo and emissivity of the surface 499 567 ! ---------------------------------------- 500 CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref) 568 DO islope = 1,nslope 569 piceco2_tmp(:) = piceco2(:,islope) 570 alb_tmp(:,:) = psolaralb(:,:,islope) 571 emisref_tmp(:) = 0. 572 CALL albedocaps(zls,ngrid,piceco2_tmp,alb_tmp,emisref_tmp) 573 psolaralb(:,1,islope) = alb_tmp(:,1) 574 psolaralb(:,2,islope) = alb_tmp(:,2) 575 emisref(:,islope) = emisref_tmp(:) 576 ENDDO 501 577 502 578 ! set pemisurf() to emissiv when there is bare surface (needed for co2snow) 503 579 DO ig=1,ngrid 504 if (piceco2(ig).eq.0) then 505 pemisurf(ig)=emissiv 506 endif 580 DO islope = 1,nslope 581 if (piceco2(ig,islope).eq.0) then 582 pemisurf(ig,islope)=emissiv 583 endif 584 ENDDO 507 585 ENDDO 508 586 … … 515 593 516 594 DO ig=1,ngrid 517 if ( condsub(ig)) then595 if (any(condsub(ig,:))) then 518 596 do l=1,nlayer 519 597 ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep … … 527 605 c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up) 528 606 c """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 529 zmflux(1) = -zcondices (ig) - zdqssed_co2(ig)607 zmflux(1) = -zcondices_mesh_avg(ig) - zdqssed_co2(ig) 530 608 DO l=1,nlayer 531 609 zmflux(l+1) = zmflux(l) - condens_layer(ig,l) … … 702 780 c CO2 snow / clouds scheme 703 781 c *************************************************************** 704 call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev, 705 & condens_layer,zcondices,zfallice,pemisurf) 782 DO islope = 1,nslope 783 emisref_tmp(:) = emisref(:,islope) 784 condsub_tmp(:) = condsub(:,islope) 785 condens_layer_tmp(:,:) = condens_layer(:,:)* 786 & cos(pi*def_slope_mean(islope)/180.) 787 zcondices_tmp(:) = zcondices(:,islope)* 788 & cos(pi*def_slope_mean(islope)/180.) 789 zfallice_tmp(:,:) = zfallice(:,:)* 790 & cos(pi*def_slope_mean(islope)/180.) 791 pemisurf_tmp(:) = pemisurf(:,islope) 792 793 call co2snow(ngrid,nlayer,ptimestep,emisref_tmp,condsub_tmp, 794 & pplev,condens_layer_tmp,zcondices_tmp,zfallice_tmp, 795 & pemisurf_tmp) 796 pemisurf(:,islope) = pemisurf_tmp(:) 797 798 ENDDO 706 799 c *************************************************************** 707 800 c Ecriture des diagnostiques … … 739 832 & 1./(bcond-acond*log(.01*vmr_co2(ngrid,1)* 740 833 & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) 741 pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep 834 DO islope = 1,nslope 835 pdtsrfc(ngrid,islope)=(ztcondsol(ngrid)- 836 & ztsrf(ngrid,islope))/ptimestep 837 ENDDO ! islope = 1,nslope 742 838 endif 743 839 endif
Note: See TracChangeset
for help on using the changeset viewer.
