Ignore:
Timestamp:
Mar 30, 2000, 5:49:22 PM (25 years ago)
Author:
lmdzadmin
Message:

Inclut la nouvelle formulation du calcul des fluxs a l'interface de JL Dufresne
LF

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ.3.3/trunk/libf/phylmd/clmain.F

    r39 r77  
    497497      REAL local_q(klon,klev)
    498498      REAL local_ts(klon)
    499       REAL zx_alf1(klon), zx_alf2(klon) !valeur ambiante par extrapola.
    500499      REAL psref(klon) ! pression de reference pour temperature potent.
    501500      REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
     
    508507      REAL zdelz
    509508c======================================================================
     509C 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)
     514c======================================================================
    510515      REAL zcor, zdelta, zcvm5
    511516#include "YOMCST.h"
     
    516521         psref(i) = paprs(i,1) !pression de reference est celle au sol
    517522         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.0
    520          zx_alf2(i) = 1.0 - zx_alf1(i)
    521523      ENDDO
    522524      DO k = 1, klev
     
    535537     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
    536538     .                 * pplay(i,1)/(RD*t(i,1))
    537          zx_coef(i,1) = zx_coef(i,1) * dtime*RG
    538539      ENDDO
    539540c
     
    557558      ENDDO
    558559c======================================================================
     560c
     561c zx_qs = qsat en kg/kg
     562c
    559563      DO i = 1, knon
    560564         IF (thermcep) THEN
     
    579583           ENDIF
    580584         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
     586c        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
     589c
    584590      ENDDO
    585591      DO i = 1, knon
     
    613619      ENDDO
    614620      ENDDO
     621C
     622C nouvelle formulation JL Dufresne
     623C
     624C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt
     625C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
     626C
    615627      DO i = 1, knon
    616628         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))
    619629         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)))
    623631     .                /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)
    625633c
    626634         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))
    629635         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)))
    633637     .                /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
     641C === Calcul de la temperature de surface ===
     642C
     643C zx_sl = chaleur latente d'evaporation ou de sublimation
     644C
     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
     651C 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)
     659c 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)
     663c 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)
     673c== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
     674c== 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)
     677c 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
    664683c==== une fois on a zx_h_ts, on peut faire l'iteration ========
    665684      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
    671687      ENDDO
    672688      DO k = 2, klev
     
    679695c== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
    680696c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
    681       DO i = 1, knon
    682         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       ENDDO
    690697      DO k = 2, klev
    691698      DO i = 1, knon
     
    698705      ENDDO
    699706c======================================================================
    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======================================================================
     707C Calcul tendances
    710708      DO k = 1, klev
    711709      DO i = 1, knon
     
    713711         d_q(i,k) = local_q(i,k) - q(i,k)
    714712      ENDDO
    715       ENDDO
    716       DO i = 1, knon
    717          d_ts(i) = local_ts(i) - ts(i)
    718713      ENDDO
    719714c
Note: See TracChangeset for help on using the changeset viewer.