Changeset 1849 for LMDZ5/trunk/libf/phylmd/fisrtilp.F90
- Timestamp:
- Aug 27, 2013, 12:55:18 PM (11 years ago)
- File:
-
- 1 edited
-
LMDZ5/trunk/libf/phylmd/fisrtilp.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phylmd/fisrtilp.F90
r1746 r1849 8 8 frac_impa, frac_nucl, beta, & 9 9 prfl, psfl, rhcl, zqta, fraca, & 10 ztv, zpspsk, ztla, zthl, iflag_cldcon) 10 ztv, zpspsk, ztla, zthl, iflag_cldcon, & 11 iflag_ice_thermo) 11 12 12 13 ! … … 51 52 REAL zpspsk(klon,klev),ztla(klon,klev) 52 53 REAL zthl(klon,klev) 54 REAL ztfondue, qsl, qsi 53 55 54 56 logical lognormale(klon) 57 logical ice_thermo 55 58 56 59 !AA … … 77 80 INTEGER ncoreczq 78 81 INTEGER iflag_cldcon 82 INTEGER iflag_ice_thermo 79 83 PARAMETER (ninter=5) 80 84 LOGICAL evap_prec ! evaporation de la pluie … … 97 101 INTEGER i, k, n, kk 98 102 REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 99 REAL zrfl(klon), zrfln(klon), zqev, zqevt 100 REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq 103 REAL zrfl(klon), zrfln(klon), zqev, zqevt 104 REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti 105 REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq 106 REAL zoliqp(klon), zoliqi(klon) 101 107 REAL ztglace, zt(klon) 102 108 INTEGER nexpo ! exponentiel pour glace/eau 103 109 REAL zdz(klon),zrho(klon),ztot , zrhol(klon) 104 110 REAL zchau ,zfroi ,zfice(klon),zneb(klon) 111 REAL zmelt, zpluie, zice, zcondold 112 PARAMETER (ztfondue=278.15) 105 113 ! 106 114 LOGICAL appel1er … … 145 153 DATA appel1er /.TRUE./ 146 154 !ym 155 ice_thermo = iflag_ice_thermo .GE. 1 147 156 zdelq=0.0 148 157 … … 188 197 ! 189 198 ztglace = RTT - 15.0 190 nexpo = 6 199 !AJ< 200 IF (ice_thermo) THEN 201 nexpo = 2 202 ELSE 203 nexpo = 6 204 ENDIF 205 !! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg) 206 !! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg) 207 !>AJ 191 208 !cc nexpo = 1 192 209 ! … … 223 240 ! DO i = 1, klon 224 241 zrfl(i) = 0.0 242 zifl(i) = 0.0 225 243 zneb(i) = seuil_neb 226 244 ENDDO … … 272 290 273 291 292 ! Calculer l'evaporation de la precipitation 293 ! 294 295 274 296 IF (evap_prec) THEN 275 297 DO i = 1, klon 276 IF (zrfl(i) .GT.0.) THEN 298 !AJ< 299 !! IF (zrfl(i) .GT.0.) THEN 300 IF (zrfl(i)+zifl(i).GT.0.) THEN 301 !>AJ 277 302 IF (thermcep) THEN 278 303 zdelta=MAX(0.,SIGN(1.,RTT-zt(i))) … … 288 313 ENDIF 289 314 ENDIF 290 zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) ) 291 zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) & 292 * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 293 zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & 294 * RG*dtime/(paprs(i,k)-paprs(i,k+1)) 295 zqev = MIN (zqev, zqevt) 296 zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & 297 /RG/dtime 298 299 ! pour la glace, on ré-évapore toute la précip dans la 300 ! couche du dessous 301 ! la glace venant de la couche du dessus est simplement 302 ! dans la couche du dessous. 303 304 IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0. 305 306 zq(i) = zq(i) - (zrfln(i)-zrfl(i)) & 307 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime 308 zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & 309 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & 310 * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) 311 zrfl(i) = zrfln(i) 312 ENDIF 313 ENDDO 314 ENDIF 315 ENDIF ! (zrfl(i)+zifl(i).GT.0.) 316 ENDDO 317 !AJ< 318 IF (.NOT. ice_thermo) THEN 319 DO i = 1, klon 320 !AJ< 321 !! IF (zrfl(i) .GT.0.) THEN 322 IF (zrfl(i)+zifl(i).GT.0.) THEN 323 !>AJ 324 zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) ) 325 zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) & 326 * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 327 zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & 328 * RG*dtime/(paprs(i,k)-paprs(i,k+1)) 329 zqev = MIN (zqev, zqevt) 330 zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & 331 /RG/dtime 332 333 ! pour la glace, on ré-évapore toute la précip dans la 334 ! couche du dessous 335 ! la glace venant de la couche du dessus est simplement 336 ! dans la couche du dessous. 337 338 IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0. 339 340 zq(i) = zq(i) - (zrfln(i)-zrfl(i)) & 341 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime 342 zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & 343 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & 344 * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) 345 zrfl(i) = zrfln(i) 346 zifl(i) = 0. 347 ENDIF ! (zrfl(i)+zifl(i).GT.0.) 348 ENDDO 349 ! 350 ELSE ! (.NOT. ice_thermo) 351 ! 352 DO i = 1, klon 353 !AJ< 354 !! IF (zrfl(i) .GT.0.) THEN 355 IF (zrfl(i)+zifl(i).GT.0.) THEN 356 !>AJ 357 !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 358 ! Modification de l'évaporation avec la glace 359 ! Différentiation entre précipitation liquide et solide 360 ! On suppose que coef_evai=2*coef_eva 361 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 362 363 zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) ) 364 ! zqev0 = MAX (0.0, zqs(i)-zq(i) ) 365 366 !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 367 ! On différencie qsat pour l'eau et la glace 368 ! Si zdelta=1. --> glace 369 ! Si zdelta=0. --> eau liquide 370 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 371 372 qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k) 373 qsl= MIN(0.5,qsl) 374 zcor= 1./(1.-RETV*qsl) 375 qsl= qsl*zcor 376 377 zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) & 378 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 379 zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & 380 *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 381 382 qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k) 383 qsi= MIN(0.5,qsi) 384 zcor= 1./(1.-RETV*qsi) 385 qsi= qsi*zcor 386 387 zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) & 388 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 389 zqevti = MAX(0.0,MIN(zqevti,zifl(i))) & 390 *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 391 392 393 !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 394 ! Vérification sur l'évaporation 395 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 396 397 IF (zqevt+zqevti.GT.zqev0) THEN 398 zqev=zqev0*zqevt/(zqevt+zqevti) 399 zqevi=zqev0*zqevti/(zqevt+zqevti) 400 401 ELSE 402 IF (zqevt+zqevti.GT.0.) THEN 403 zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt) 404 zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti) 405 ELSE 406 zqev=0. 407 zqevi=0. 408 ENDIF 409 ENDIF 410 411 zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & 412 /RG/dtime) 413 zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) & 414 /RG/dtime) 415 416 ! Pour la glace, on révapore toute la précip dans la couche du dessous 417 ! la glace venant de la couche du dessus est simplement dans la couche 418 ! du dessous. 419 420 ! IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0. 421 ! print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip' 422 zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & 423 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime 424 zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & 425 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & 426 * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) & 427 + (zifln(i)-zifl(i)) & 428 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & 429 * RLSTT/RCPD/(1.0+RVTMP2*zq(i)) 430 431 zrfl(i) = zrfln(i) 432 zifl(i) = zifln(i) 433 434 ENDIF ! (zrfl(i)+zifl(i).GT.0.) 435 ENDDO 436 437 ENDIF ! (.NOT. ice_thermo) 438 439 ENDIF ! (evap_prec) 315 440 ! 316 441 ! Calculer Qs et L/Cp*dQs/dT: … … 478 603 zq(i) = zq(i) - zcond(i) 479 604 ! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD 480 zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) 481 ENDDO 605 ENDDO 606 !AJ< 607 IF (.NOT. ice_thermo) THEN 608 DO i = 1, klon 609 zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) 610 ENDDO 611 ELSE 612 DO i = 1, klon 613 zfice(i) = 1.0 - (zt(i)-ztglace) / (273.15-ztglace) 614 zfice(i) = MIN(MAX(zfice(i),0.0),1.0) 615 zfice(i) = zfice(i)**nexpo 616 zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) & 617 +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i)) 618 ! print*,zt(i),zrfl(i),zifl(i),'temp1' 619 ENDDO 620 ENDIF 621 !>AJ 482 622 ! 483 623 ! Partager l'eau condensee en precipitation et eau liquide nuageuse … … 488 628 zrho(i) = pplay(i,k) / zt(i) / RD 489 629 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG) 630 ENDIF 631 ENDDO 632 !AJ< 633 IF (.NOT. ice_thermo) THEN 634 DO i = 1, klon 635 IF (rneb(i,k).GT.0.0) THEN 490 636 zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace) 491 637 zfice(i) = MIN(MAX(zfice(i),0.0),1.0) 492 638 zfice(i) = zfice(i)**nexpo 639 !! zfice(i)=0. 640 ENDIF 641 ENDDO 642 ENDIF 643 DO i = 1, klon 644 IF (rneb(i,k).GT.0.0) THEN 493 645 zneb(i) = MAX(rneb(i,k), seuil_neb) 646 ! zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) 647 ! print*,zt(i),'fractionglace' 648 !>AJ 494 649 radliq(i,k) = zoliq(i)/REAL(ninter+1) 495 650 ENDIF … … 502 657 503 658 IF (zneb(i).EQ.seuil_neb) THEN 659 zpluie= 0.0 660 zice = 0.0 504 661 ztot = 0.0 505 662 ELSE … … 519 676 zchau = zct *dtime/REAL(ninter) * zoliq(i) & 520 677 *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )**2)) *(1.-zfice(i)) 521 ztot = zchau + zfroi 678 !AJ< 679 IF (.NOT. ice_thermo) THEN 680 ztot = zchau + zfroi 681 ELSE 682 zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i))) 683 zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) 684 ztot = zpluie + zice 685 ENDIF 686 !>AJ 522 687 ztot = MAX(ztot ,0.0) 523 688 ENDIF 524 689 ztot = MIN(ztot,zoliq(i)) 690 !AJ< 691 ! zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0) 692 ! zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice , 0.0) 693 zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0) 694 zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice , 0.0) 525 695 zoliq(i) = MAX(zoliq(i)-ztot , 0.0) 696 !>AJ 526 697 radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1) 527 698 ENDIF … … 529 700 ENDDO 530 701 ! 531 DO i = 1, klon 532 IF (rneb(i,k).GT.0.0) THEN 702 IF (.NOT. ice_thermo) THEN 703 DO i = 1, klon 704 IF (rneb(i,k).GT.0.0) THEN 533 705 d_ql(i,k) = zoliq(i) 534 706 zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) & 535 707 * (paprs(i,k)-paprs(i,k+1))/(RG*dtime) 536 ENDIF 537 IF (zt(i).LT.RTT) THEN 708 ENDIF 709 ENDDO 710 ELSE 711 DO i = 1, klon 712 IF (rneb(i,k).GT.0.0) THEN 713 d_ql(i,k) = zoliq(i) 714 !AJ< 715 zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) & 716 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 717 zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) & 718 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 719 ! zrfl(i) = zrfl(i)+ zpluie & 720 ! *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 721 ! zifl(i) = zifl(i)+ zice & 722 ! *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 723 724 ENDIF 725 ENDDO 726 ENDIF 727 728 IF (ice_thermo) THEN 729 DO i = 1, klon 730 zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2 731 zmelt = MIN(MAX(zmelt,0.),1.) 732 zrfl(i)=zrfl(i)+zmelt*zifl(i) 733 zifl(i)=zifl(i)*(1.-zmelt) 734 ! print*,zt(i),'octavio1' 735 zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) & 736 *RLMLT/RCPD/(1.0+RVTMP2*zq(i)) 737 ! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2' 738 ENDDO 739 ENDIF 740 741 742 IF (.NOT. ice_thermo) THEN 743 DO i = 1, klon 744 IF (zt(i).LT.RTT) THEN 538 745 psfl(i,k)=zrfl(i) 539 ELSE746 ELSE 540 747 prfl(i,k)=zrfl(i) 541 ENDIF 542 ENDDO 748 ENDIF 749 ENDDO 750 ELSE 751 ! JAM************************************************* 752 ! Revoir partie ci-dessous: à quoi servent psfl et prfl? 753 ! ***************************************************** 754 755 DO i = 1, klon 756 ! IF (zt(i).LT.RTT) THEN 757 psfl(i,k)=zifl(i) 758 ! ELSE 759 prfl(i,k)=zrfl(i) 760 ! ENDIF 761 !>AJ 762 ENDDO 763 ENDIF 764 ! 543 765 ! 544 766 ! Calculer les tendances de q et de t: … … 604 826 DO i = 1, klon 605 827 IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN 606 snow(i) = zrfl(i) 828 !AJ< 829 !! snow(i) = zrfl(i) 830 snow(i) = zrfl(i)+zifl(i) 831 !>AJ 607 832 zlh_solid(i) = RLSTT-RLVTT 608 833 ELSE
Note: See TracChangeset
for help on using the changeset viewer.
