Changeset 77 for LMDZ.3.3/trunk/libf
- Timestamp:
- Mar 30, 2000, 5:49:22 PM (25 years ago)
- File:
-
- 1 edited
-
LMDZ.3.3/trunk/libf/phylmd/clmain.F (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ.3.3/trunk/libf/phylmd/clmain.F
r39 r77 497 497 REAL local_q(klon,klev) 498 498 REAL local_ts(klon) 499 REAL zx_alf1(klon), zx_alf2(klon) !valeur ambiante par extrapola.500 499 REAL psref(klon) ! pression de reference pour temperature potent. 501 500 REAL zx_pkh(klon,klev), zx_pkf(klon,klev) … … 508 507 REAL zdelz 509 508 c====================================================================== 509 C Variables intermediaires pour le calcul des fluxs a la surface 510 real zx_mh(klon), zx_nh(klon), zx_oh(klon) 511 real zx_mq(klon), zx_nq(klon), zx_oq(klon) 512 real zx_k1(klon), zx_dq_s_dt(klon) 513 real zx_qsat(klon) 514 c====================================================================== 510 515 REAL zcor, zdelta, zcvm5 511 516 #include "YOMCST.h" … … 516 521 psref(i) = paprs(i,1) !pression de reference est celle au sol 517 522 local_ts(i) = ts(i) 518 ccc zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))519 zx_alf1(i) = 1.0520 zx_alf2(i) = 1.0 - zx_alf1(i)521 523 ENDDO 522 524 DO k = 1, klev … … 535 537 . * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) 536 538 . * pplay(i,1)/(RD*t(i,1)) 537 zx_coef(i,1) = zx_coef(i,1) * dtime*RG538 539 ENDDO 539 540 c … … 557 558 ENDDO 558 559 c====================================================================== 560 c 561 c zx_qs = qsat en kg/kg 562 c 559 563 DO i = 1, knon 560 564 IF (thermcep) THEN … … 579 583 ENDIF 580 584 ENDIF 581 c 582 zx_dq0(i) = zx_dq_s_dh 583 zx_cq0(i) = zx_qs -zx_dq_s_dh *local_ts(i)*RCPD*zx_pkh(i,1) 585 586 c zx_dq_s_dh = 0. 587 zx_dq_s_dt(i) = RCPD * zx_pkh(i,1) * zx_dq_s_dh 588 zx_qsat(i) = zx_qs 589 c 584 590 ENDDO 585 591 DO i = 1, knon … … 613 619 ENDDO 614 620 ENDDO 621 C 622 C nouvelle formulation JL Dufresne 623 C 624 C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt 625 C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt 626 C 615 627 DO i = 1, knon 616 628 zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2)) 617 . + zx_coef(i,1)*beta(i)618 . *(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))619 629 zx_cq(i,1) = (local_q(i,1)*delp(i,1) 620 . +zx_coef(i,2)*z_gamaq(i,2) 621 . +zx_cq(i,2)*(zx_coef(i,2) 622 . -zx_coef(i,1)*zx_alf2(i)*beta(i)) ) 630 . +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2))) 623 631 . /zx_buf1(i) 624 zx_dq(i,1) = beta(i) * zx_coef(i,1)/ zx_buf1(i)632 zx_dq(i,1) = -1. * RG / zx_buf1(i) 625 633 c 626 634 zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2)) 627 . +zx_coef(i,1)628 . *(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))629 635 zx_ch(i,1) = (local_h(i,1)*delp(i,1) 630 . +zx_coef(i,2)*z_gamah(i,2) 631 . +zx_ch(i,2)*(zx_coef(i,2) 632 . -zx_coef(i,1)*zx_alf2(i)) ) 636 . +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2))) 633 637 . /zx_buf2(i) 634 zx_dh(i,1) = zx_coef(i,1) / zx_buf2(i) 635 ENDDO 636 c==== calculer d'abord local_h(0) ============================ 637 DO i = 1, knon 638 zx_sl(i) = RLVTT 639 if (local_ts(i) .LT. RTT) zx_sl(i) = RLSTT 640 ENDDO 641 DO i = 1, knon 642 zx_h_grnd = t_grnd * zx_pkh(i,1) *RCPD 643 zx_h_ts(i) = local_ts(i) * zx_pkh(i,1) *RCPD 644 zx_a = radsol(i) 645 . + zx_coef(i,1)/(RG*dtime)*beta(i)*zx_sl(i) 646 . *(zx_cq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2)) 647 . +zx_alf2(i)*zx_cq(i,2)) 648 . + zx_coef(i,1)/(RG*dtime)*beta(i)*zx_sl(i) 649 . *((zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2)) 650 . -1.0)*zx_cq0(i)) 651 . + zx_coef(i,1)/(RG*dtime)/zx_pkh(i,1) 652 . *(zx_ch(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2)) 653 . +zx_alf2(i)*zx_ch(i,2)) 654 zx_a = zx_h_ts(i) + dif_grnd(i)*dtime*zx_h_grnd 655 . + cal(i)*zx_pkh(i,1)*dtime * zx_a 656 zx_b = zx_coef(i,1)/(RG*dtime) *beta(i)*zx_sl(i)*zx_dq0(i) 657 . *(1.0-zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))) 658 . + zx_coef(i,1)/(RG*dtime)/zx_pkh(i,1) 659 . *(1.0-zx_dh(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))) 660 zx_b = 1.0 + dif_grnd(i)*dtime + cal(i)*zx_pkh(i,1)*dtime*zx_b 661 zx_h_ts(i) = zx_a / zx_b 662 local_ts(i) = zx_h_ts(i) / zx_pkh(i,1)/RCPD 663 ENDDO 638 zx_dh(i,1) = -1. * RG / zx_buf2(i) 639 ENDDO 640 641 C === Calcul de la temperature de surface === 642 C 643 C zx_sl = chaleur latente d'evaporation ou de sublimation 644 C 645 do i = 1, knon 646 zx_sl(i) = RLVTT 647 if (local_ts(i) .LT. RTT) zx_sl(i) = RLSTT 648 zx_k1(i) = zx_coef(i,1) 649 enddo 650 do i = 1, knon 651 C Q 652 zx_oq(i) = 1. - (beta(i) * zx_k1(i) * zx_dq(i,1) * dtime) 653 zx_mq(i) = beta(i) * zx_k1(i) * 654 . (zx_cq(i,1) - zx_qsat(i) 655 . + zx_dq_s_dt(i) * local_ts(i)) 656 . / zx_oq(i) 657 zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) 658 . / zx_oq(i) 659 c H 660 zx_oh(i) = 1. - (zx_k1(i) * zx_dh(i,1) * dtime) 661 zx_mh(i) = zx_k1(i) * zx_ch(i,1) / zx_oh(i) 662 zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i,1))/ zx_oh(i) 663 c Tsurface 664 local_ts(i) = (ts(i) + cal(i)/(RCPD * zx_pkh(i,1)) * dtime * 665 . (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) 666 . + dif_grnd(i) * t_grnd * dtime)/ 667 . ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i,1)) * ( 668 . zx_nh(i) + zx_sl(i) * zx_nq(i)) 669 . + dtime * dif_grnd(i)) 670 zx_h_ts(i) = local_ts(i) * RCPD * zx_pkh(i,1) 671 d_ts(i) = local_ts(i) - ts(i) 672 zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i) 673 c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas 674 c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) 675 flux_q(i,1) = zx_mq(i) + zx_nq(i) * local_ts(i) 676 flux_t(i,1) = zx_mh(i) + zx_nh(i) * local_ts(i) 677 c Derives des flux dF/dTs (W m-2 K-1): 678 dflux_s(i) = zx_nh(i) 679 dflux_l(i) = (zx_sl(i) * zx_nq(i)) 680 ENDDO 681 682 664 683 c==== une fois on a zx_h_ts, on peut faire l'iteration ======== 665 684 DO i = 1, knon 666 zx_q_0(i) = zx_cq0(i) + zx_dq0(i)*zx_h_ts(i) 667 ENDDO 668 DO i = 1, knon 669 local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*zx_h_ts(i) 670 local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*zx_q_0(i) 685 local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i)*dtime 686 local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*flux_q(i)*dtime 671 687 ENDDO 672 688 DO k = 2, klev … … 679 695 c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas 680 696 c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) 681 DO i = 1, knon682 flux_q(i,1) = (zx_coef(i,1)/dtime/RG) * beta(i)683 . * (local_q(i,1)*zx_alf1(i)684 . +local_q(i,2)*zx_alf2(i)-zx_q_0(i))685 flux_t(i,1) = (zx_coef(i,1)/dtime/RG)686 . * (local_h(i,1)*zx_alf1(i)687 . +local_h(i,2)*zx_alf2(i)-zx_h_ts(i))688 . / zx_pkh(i,1)689 ENDDO690 697 DO k = 2, klev 691 698 DO i = 1, knon … … 698 705 ENDDO 699 706 c====================================================================== 700 c Derives des flux dF/dTs (W m-2 K-1): 701 DO i = 1, knon 702 dflux_s(i) = (RCPD*zx_coef(i,1))/(RG*dtime)/zx_pkh(i,1) 703 . *(zx_dh(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))-1.0) 704 dflux_l(i) = (beta(i)*RLVTT*RCPD*zx_coef(i,1))/(RG*dtime) 705 . *zx_dq0(i) 706 . *(zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))-1.0) 707 dflux_g(i) = dif_grnd(i) * RCPD 708 ENDDO 709 c====================================================================== 707 C Calcul tendances 710 708 DO k = 1, klev 711 709 DO i = 1, knon … … 713 711 d_q(i,k) = local_q(i,k) - q(i,k) 714 712 ENDDO 715 ENDDO716 DO i = 1, knon717 d_ts(i) = local_ts(i) - ts(i)718 713 ENDDO 719 714 c
Note: See TracChangeset
for help on using the changeset viewer.
