Changeset 473
- Timestamp:
- Dec 14, 2011, 1:20:09 PM (13 years ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 2 edited
-
README (modified) (1 diff)
-
libf/phymars/vdifc.F (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r472 r473 1339 1339 order of tracers was changed). 1340 1340 1341 == 14/12/2011 == AC 1342 >> Improved computation of mixing in vdifc for temperatures less than the saturation 1343 temperature. 1344 The scheme works by computing mixing of T and then liberating latent heat when T<Ts. 1345 Taking into account the mass of co2 ice condensated, a new mixing profile for T is computed. 1346 When the computation has converged, the tendancy for mixing of T is taken with respect to 1347 the initial T profile for which co2 ice has been condensated. 1348 >> This scheme is based on initial work done by Melanie Vangvichith and has been adapted/tested/bug-fixed 1349 and improved to Mars. 1350 -
trunk/LMDZ.MARS/libf/phymars/vdifc.F
r325 r473 72 72 c ------ 73 73 74 REAL ust(ngridmx),tst(ngridmx) 75 74 76 INTEGER ilev,ig,ilay,nlev 75 77 … … 112 114 REAL kmixmin 113 115 116 117 c Mass-variation scheme : 118 c ~~~~~~~ 119 120 INTEGER j,l 121 REAL zcondicea(ngrid,nlayermx) 122 REAL zt(ngridmx,nlayermx),ztcond(ngridmx,nlayermx+1) 123 REAL betam(ngridmx,nlayermx),dmice(ngridmx,nlayermx) 124 REAL pdtc(ngrid,nlayermx) 125 REAL zhs(ngridmx,nlayermx) 126 REAL cpice,ccond 127 SAVE ccond,cpice 128 DATA cpice /1000./ 129 130 c Theta_m formulation for mass-variation scheme : 131 c ~~~~~~~ 132 133 INTEGER ico2,llnt(ngridmx) 134 SAVE ico2 135 REAL m_co2, m_noco2, A , B 136 SAVE A, B, m_co2, m_noco2 137 REAL vmr_co2(ngridmx,nlayermx) 138 REAL qco2,mmean 139 114 140 c ** un petit test de coherence 115 141 c -------------------------- … … 126 152 bcond=1./tcond1mb 127 153 acond=r/latcond 154 ccond=cpp/(g*latcond) 128 155 PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond 129 PRINT*,' acond,bcond',acond,bcond 156 PRINT*,' acond,bcond,ccond',acond,bcond,ccond 157 158 159 ico2=0 160 161 if (tracer) then 162 c Prepare Special treatment if one of the tracer is CO2 gas 163 do iq=1,nqmx 164 if (noms(iq).eq."co2") then 165 ico2=iq 166 m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) 167 m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) 168 c Compute A and B coefficient use to compute 169 c mean molecular mass Mair defined by 170 c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2 171 c 1/Mair = A*q(ico2) + B 172 A =(1/m_co2 - 1/m_noco2) 173 B=1/m_noco2 174 endif 175 enddo 176 end if 130 177 131 178 firstcall=.false. … … 149 196 DO ig=1,ngrid 150 197 za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g 198 ! Mass variation scheme: 199 betam(ig,ilay)=-za(ig,ilay)*latcond/(cpp*ppopsk(ig,ilay)) 151 200 ENDDO 152 201 ENDDO … … 163 212 ENDDO 164 213 DO ig=1,ngrid 165 214 zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig)) 166 215 ENDDO 167 216 … … 185 234 ENDIF 186 235 236 c ----------------------------------- 187 237 c Potential Condensation temperature: 188 238 c ----------------------------------- 189 239 190 c if (callcond) then 191 c DO ilev=1,nlay 192 c DO ig=1,ngrid 193 c zhcond(ig,ilev) = 194 c & (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev) 195 c END DO 196 c END DO 197 c else 198 call zerophys(ngrid*nlay,zhcond) 199 c end if 240 c Compute CO2 Volume mixing ratio 241 c ------------------------------- 242 if (callcond.and.(ico2.ne.0)) then 243 DO ilev=1,nlay 244 DO ig=1,ngrid 245 qco2=pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep 246 c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2) 247 mmean=1/(A*qco2 +B) 248 vmr_co2(ig,ilev) = qco2*mmean/m_co2 249 ENDDO 250 ENDDO 251 else 252 DO ilev=1,nlay 253 DO ig=1,ngrid 254 vmr_co2(ig,ilev)=0.95 255 ENDDO 256 ENDDO 257 end if 258 259 c forecast of atmospheric temperature zt and frost temperature ztcond 260 c -------------------------------------------------------------------- 261 262 DO ilev=1,nlay 263 DO ig=1,ngrid 264 ztcond(ig,ilev)= 265 & 1./(bcond-acond*log(.01*vmr_co2(ig,ilev)*pplay(ig,ilev))) 266 if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica 267 ENDDO 268 ENDDO 269 270 ztcond(:,nlay+1)=ztcond(:,nlay) 271 272 if (callcond) then 273 DO ilev=1,nlay 274 DO ig=1,ngrid 275 ! zhcond(ig,ilev) = 276 ! & (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev) 277 zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev) 278 END DO 279 END DO 280 else 281 call zerophys(ngrid*nlay,zhcond) 282 end if 200 283 201 284 … … 209 292 zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep 210 293 zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep 211 zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))294 ! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev)) 212 295 ENDDO 213 296 ENDDO … … 239 322 zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+(0.7*wmax(:))**2) ! subgrid gustiness added by quadratic interpolation with a coeff beta 240 323 zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2) ! LES comparisons. This parameter is linked to the thermals settings) 324 ! ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2) 325 ! tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:) 326 ! ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)) 327 ! tst(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:)) 241 328 ELSE 242 329 zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance 243 330 zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance 331 ! ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)) 332 ! tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:) 244 333 ENDIF 245 334 … … 254 343 ! & 2,zcdh_true) 255 344 ! call WRITEDIAGFI(ngridmx,'u*', 256 ! & 'friction velocity','m/s', 257 ! & 2,sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2)) 258 ! call WRITEDIAGFI(ngridmx,'T*', 259 ! & 'friction temperature','K', 260 ! & 2,zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2)* 261 ! & -(ptsrf(:)-ph(:,1)) 262 ! & /(sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2))) 345 ! & 'friction velocity','m/s',2,ust) 346 ! call WRITEDIAGFI(ngridmx,'T*', 347 ! & 'friction temperature','K',2,tst) 263 348 ! call WRITEDIAGFI(ngridmx,'rm-1', 264 349 ! & 'aerodyn momentum conductance','m/s', … … 405 490 c ------------- 406 491 492 c Mass variation scheme: 407 493 CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2)) 408 494 CALL multipl(ngrid,zcdh,zb0,zb) 409 495 410 DO ig=1,ngrid 411 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 412 zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig) 413 zd(ig,nlay)=zb(ig,nlay)*z1(ig) 414 ENDDO 415 416 DO ilay=nlay-1,1,-1 417 DO ig=1,ngrid 418 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 419 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 420 zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ 421 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) 422 zd(ig,ilay)=zb(ig,ilay)*z1(ig) 423 ENDDO 424 ENDDO 496 c on initialise dm c 497 498 pdtc(:,:)=0. 499 zt(:,:)=0. 500 dmice(:,:)=0. 425 501 426 502 c ** calcul de (d Planck / dT) a la temperature d'interface … … 431 507 zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig) 432 508 ENDDO 509 510 511 ! calcul de zc et zd pour la couche top en prenant en compte le terme 512 ! de variation de masse (on fait une boucle pour que ça converge) 513 514 ! Identification des points de grilles qui ont besoin de la correction 515 516 llnt(:)=1 517 518 DO ig=1,ngrid 519 DO l=1,nlay 520 if(zh(ig,l) .lt. zhcond(ig,l)) then 521 llnt(ig)=300 522 ! 200 and 100 do not go beyond month 9 with normal dissipation 523 goto 5 524 endif 525 ENDDO 526 5 continue 527 ENDDO 528 529 DO ig=1,ngrid 530 531 ! Initialization of z1 and zd, which do not depend on dmice 532 533 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 534 zd(ig,nlay)=zb(ig,nlay)*z1(ig) 535 536 DO ilay=nlay-1,1,-1 537 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 538 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 539 zd(ig,ilay)=zb(ig,ilay)*z1(ig) 540 ENDDO 541 542 ! Convergence loop 543 544 DO j=1,llnt(ig) 545 546 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 547 zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay) 548 & -betam(ig,nlay)*dmice(ig,nlay) 549 zc(ig,nlay)=zc(ig,nlay)*z1(ig) 550 ! zd(ig,nlay)=zb(ig,nlay)*z1(ig) 551 552 ! calcul de zc et zd pour les couches du haut vers le bas 553 554 DO ilay=nlay-1,1,-1 555 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 556 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 557 zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ 558 $ zb(ig,ilay+1)*zc(ig,ilay+1)- 559 $ betam(ig,ilay)*dmice(ig,ilay))*z1(ig) 560 ! zd(ig,ilay)=zb(ig,ilay)*z1(ig) 561 ENDDO 433 562 434 563 c ** calcul de la temperature_d'interface et de sa tendance. … … 444 573 c ---------------------------------------------------- 445 574 446 DO ig=1,ngrid447 575 z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) 448 576 s +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep 449 577 z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig) 450 578 ztsrf2(ig)=z1(ig)/z2(ig) 451 pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep 452 453 c Modif speciale CO2 condensation: 454 c tconds = 1./(bcond-acond*log(.0095*pplev(ig,1))) 455 c if ((callcond).and. 456 c & ((co2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then 457 c zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds 458 c else 459 zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) 460 c end if 461 ENDDO 579 ! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop 580 581 zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) 462 582 463 583 c ** et a partir de la temperature au sol on remonte 464 584 c ----------------------------------------------- 465 585 466 DO ilay=2,nlay 467 DO ig=1,ngrid 468 hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif co2cond 469 zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh 470 ENDDO 471 ENDDO 586 DO ilay=2,nlay 587 zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1) 588 ENDDO 589 DO ilay=1,nlay 590 zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay) 591 ENDDO 592 593 c Condensation/sublimation in the atmosphere 594 c ------------------------------------------ 595 c (computation of zcondicea and dmice) 596 597 zcondicea(ig,:)=0. 598 pdtc(ig,:)=0. 599 600 DO l=nlay , 1, -1 601 IF(zt(ig,l).LT.ztcond(ig,l)) THEN 602 pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep 603 zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1)) 604 & *ccond*pdtc(ig,l) 605 dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep 606 END IF 607 ENDDO 608 609 ENDDO !of Do j=1,XXX 610 611 ENDDO !of Do ig=1,nlay 612 613 pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep 472 614 473 615 … … 570 712 ! special case for water ice tracer: do not include 571 713 ! h2o ice tracer from surface (which is set when handling 572 714 ! h2o vapour case (see further down). 573 715 DO ig=1,ngrid 574 716 z1(ig)=1./(za(ig,1)+zb(ig,1)+ … … 641 783 pdvdif(ig,ilev)=( zv(ig,ilev)- 642 784 $ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep 643 hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep , 644 $ zhcond(ig,ilev)) ! modif co2cond 645 pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep 646 ENDDO 647 ENDDO 648 649 785 hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep 786 $ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev) 787 pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep 788 ENDDO 789 ENDDO 790 650 791 if (tracer) then 651 792 DO iq = 1, nq … … 670 811 DO ilev=1,nlay 671 812 WRITE(*,'(4f15.7)') 672 s ph(ngrid/2+1,ilev),zh (ngrid/2+1,ilev),813 s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev), 673 814 s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev) 674 815 675 816 ENDDO 676 817 ENDIF 677 678 if (calltherm .and. outptherm) then679 if (ngrid .eq. 1) then680 call WRITEDIAGFI(ngrid,'zh','zh inside vdifc',681 & 'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep)682 call WRITEDIAGFI(ngrid,'zkh','zkh',683 & 'SI',1,zkh)684 endif685 endif686 818 687 819 RETURN
Note: See TracChangeset
for help on using the changeset viewer.
