Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F (revision 83) +++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F (revision 84) @@ -497,5 +497,4 @@ REAL local_q(klon,klev) REAL local_ts(klon) - REAL zx_alf1(klon), zx_alf2(klon) !valeur ambiante par extrapola. REAL psref(klon) ! pression de reference pour temperature potent. REAL zx_pkh(klon,klev), zx_pkf(klon,klev) @@ -508,4 +507,10 @@ REAL zdelz c====================================================================== +C Variables intermediaires pour le calcul des fluxs a la surface + real zx_mh(klon), zx_nh(klon), zx_oh(klon) + real zx_mq(klon), zx_nq(klon), zx_oq(klon) + real zx_k1(klon), zx_dq_s_dt(klon) + real zx_qsat(klon) +c====================================================================== REAL zcor, zdelta, zcvm5 #include "YOMCST.h" @@ -516,7 +521,4 @@ psref(i) = paprs(i,1) !pression de reference est celle au sol local_ts(i) = ts(i) -ccc zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2)) - zx_alf1(i) = 1.0 - zx_alf2(i) = 1.0 - zx_alf1(i) ENDDO DO k = 1, klev @@ -535,5 +537,4 @@ . * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) . * pplay(i,1)/(RD*t(i,1)) - zx_coef(i,1) = zx_coef(i,1) * dtime*RG ENDDO c @@ -557,4 +558,7 @@ ENDDO c====================================================================== +c +c zx_qs = qsat en kg/kg +c DO i = 1, knon IF (thermcep) THEN @@ -579,7 +583,9 @@ ENDIF ENDIF -c - zx_dq0(i) = zx_dq_s_dh - zx_cq0(i) = zx_qs -zx_dq_s_dh *local_ts(i)*RCPD*zx_pkh(i,1) + +c zx_dq_s_dh = 0. + zx_dq_s_dt(i) = RCPD * zx_pkh(i,1) * zx_dq_s_dh + zx_qsat(i) = zx_qs +c ENDDO DO i = 1, knon @@ -613,60 +619,70 @@ ENDDO ENDDO +C +C nouvelle formulation JL Dufresne +C +C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt +C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt +C DO i = 1, knon zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2)) - . + zx_coef(i,1)*beta(i) - . *(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2)) zx_cq(i,1) = (local_q(i,1)*delp(i,1) - . +zx_coef(i,2)*z_gamaq(i,2) - . +zx_cq(i,2)*(zx_coef(i,2) - . -zx_coef(i,1)*zx_alf2(i)*beta(i)) ) + . +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2))) . /zx_buf1(i) - zx_dq(i,1) = beta(i) * zx_coef(i,1) / zx_buf1(i) + zx_dq(i,1) = -1. * RG / zx_buf1(i) c zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2)) - . +zx_coef(i,1) - . *(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2)) zx_ch(i,1) = (local_h(i,1)*delp(i,1) - . +zx_coef(i,2)*z_gamah(i,2) - . +zx_ch(i,2)*(zx_coef(i,2) - . -zx_coef(i,1)*zx_alf2(i)) ) + . +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2))) . /zx_buf2(i) - zx_dh(i,1) = zx_coef(i,1) / zx_buf2(i) - ENDDO -c==== calculer d'abord local_h(0) ============================ - DO i = 1, knon - zx_sl(i) = RLVTT - if (local_ts(i) .LT. RTT) zx_sl(i) = RLSTT - ENDDO - DO i = 1, knon - zx_h_grnd = t_grnd * zx_pkh(i,1) *RCPD - zx_h_ts(i) = local_ts(i) * zx_pkh(i,1) *RCPD - zx_a = radsol(i) - . + zx_coef(i,1)/(RG*dtime)*beta(i)*zx_sl(i) - . *(zx_cq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2)) - . +zx_alf2(i)*zx_cq(i,2)) - . + zx_coef(i,1)/(RG*dtime)*beta(i)*zx_sl(i) - . *((zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2)) - . -1.0)*zx_cq0(i)) - . + zx_coef(i,1)/(RG*dtime)/zx_pkh(i,1) - . *(zx_ch(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2)) - . +zx_alf2(i)*zx_ch(i,2)) - zx_a = zx_h_ts(i) + dif_grnd(i)*dtime*zx_h_grnd - . + cal(i)*zx_pkh(i,1)*dtime * zx_a - zx_b = zx_coef(i,1)/(RG*dtime) *beta(i)*zx_sl(i)*zx_dq0(i) - . *(1.0-zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))) - . + zx_coef(i,1)/(RG*dtime)/zx_pkh(i,1) - . *(1.0-zx_dh(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))) - zx_b = 1.0 + dif_grnd(i)*dtime + cal(i)*zx_pkh(i,1)*dtime*zx_b - zx_h_ts(i) = zx_a / zx_b - local_ts(i) = zx_h_ts(i) / zx_pkh(i,1)/RCPD - ENDDO + zx_dh(i,1) = -1. * RG / zx_buf2(i) + ENDDO + +C === Calcul de la temperature de surface === +C +C zx_sl = chaleur latente d'evaporation ou de sublimation +C + do i = 1, knon + zx_sl(i) = RLVTT + if (local_ts(i) .LT. RTT) zx_sl(i) = RLSTT + zx_k1(i) = zx_coef(i,1) + enddo + do i = 1, knon +C Q + zx_oq(i) = 1. - (beta(i) * zx_k1(i) * zx_dq(i,1) * dtime) + zx_mq(i) = beta(i) * zx_k1(i) * + . (zx_cq(i,1) - zx_qsat(i) + . + zx_dq_s_dt(i) * local_ts(i)) + . / zx_oq(i) + zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) + . / zx_oq(i) +c H + zx_oh(i) = 1. - (zx_k1(i) * zx_dh(i,1) * dtime) + zx_mh(i) = zx_k1(i) * zx_ch(i,1) / zx_oh(i) + zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i,1))/ zx_oh(i) +c Tsurface + local_ts(i) = (ts(i) + cal(i)/(RCPD * zx_pkh(i,1)) * dtime * + . (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) + . + dif_grnd(i) * t_grnd * dtime)/ + . ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i,1)) * ( + . zx_nh(i) + zx_sl(i) * zx_nq(i)) + . + dtime * dif_grnd(i)) + zx_h_ts(i) = local_ts(i) * RCPD * zx_pkh(i,1) + d_ts(i) = local_ts(i) - ts(i) + zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i) +c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas +c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) + flux_q(i,1) = zx_mq(i) + zx_nq(i) * local_ts(i) + flux_t(i,1) = zx_mh(i) + zx_nh(i) * local_ts(i) +c Derives des flux dF/dTs (W m-2 K-1): + dflux_s(i) = zx_nh(i) + dflux_l(i) = (zx_sl(i) * zx_nq(i)) + ENDDO + + c==== une fois on a zx_h_ts, on peut faire l'iteration ======== DO i = 1, knon - zx_q_0(i) = zx_cq0(i) + zx_dq0(i)*zx_h_ts(i) - ENDDO - DO i = 1, knon - local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*zx_h_ts(i) - local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*zx_q_0(i) + local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i)*dtime + local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*flux_q(i)*dtime ENDDO DO k = 2, klev @@ -679,13 +695,4 @@ c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) - DO i = 1, knon - flux_q(i,1) = (zx_coef(i,1)/dtime/RG) * beta(i) - . * (local_q(i,1)*zx_alf1(i) - . +local_q(i,2)*zx_alf2(i)-zx_q_0(i)) - flux_t(i,1) = (zx_coef(i,1)/dtime/RG) - . * (local_h(i,1)*zx_alf1(i) - . +local_h(i,2)*zx_alf2(i)-zx_h_ts(i)) - . / zx_pkh(i,1) - ENDDO DO k = 2, klev DO i = 1, knon @@ -698,14 +705,5 @@ ENDDO c====================================================================== -c Derives des flux dF/dTs (W m-2 K-1): - DO i = 1, knon - dflux_s(i) = (RCPD*zx_coef(i,1))/(RG*dtime)/zx_pkh(i,1) - . *(zx_dh(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))-1.0) - dflux_l(i) = (beta(i)*RLVTT*RCPD*zx_coef(i,1))/(RG*dtime) - . *zx_dq0(i) - . *(zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))-1.0) - dflux_g(i) = dif_grnd(i) * RCPD - ENDDO -c====================================================================== +C Calcul tendances DO k = 1, klev DO i = 1, knon @@ -713,7 +711,4 @@ d_q(i,k) = local_q(i,k) - q(i,k) ENDDO - ENDDO - DO i = 1, knon - d_ts(i) = local_ts(i) - ts(i) ENDDO c