Changeset 385 for LMDZ.3.3/branches


Ignore:
Timestamp:
Jul 15, 2002, 1:29:54 PM (22 years ago)
Author:
lmdzadmin
Message:

Definition d'un temps de relaxation pour le schema de nuages
LF

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F

    r373 r385  
    551551      save ok_newmicro
    552552      save fact_cldcon,facttemps
     553      real facteur
    553554
    554555      integer iflag_cldcon
    … …  
    578579c
    579580      REAL t_seri(klon,klev), q_seri(klon,klev)
    580       REAL ql_seri(klon,klev)
     581      REAL ql_seri(klon,klev),qs_seri(klon,klev)
    581582      REAL u_seri(klon,klev), v_seri(klon,klev)
    582583c
    … …  
    619620      character*40 vartitle
    620621      character*20 varunits
     622C     Variables liees au bilan d'energie et d'enthalpi
     623      INTEGER   if_ebil ! level for energy conserv. dignostics
     624      SAVE      if_ebil
     625      REAL ztsol(klon)
     626      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
     627     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
     628      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
     629     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
     630      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
     631      REAL      d_h_vcol_phy
     632      REAL      fs_bound, fq_bound
     633      SAVE      d_h_vcol_phy
     634      REAL      zero_v(klon)
     635      CHARACTER*15 ztit
     636      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
     637      SAVE      ip_ebil
     638      DATA      ip_ebil/2/
    621639c
    622640c Declaration des constantes et des fonctions thermodynamiques
    … …  
    627645c======================================================================
    628646      modname = 'physiq'
     647      if_ebil = 2
     648      IF (if_ebil.ge.1) THEN
     649        DO i=1,klon
     650          zero_v(i)=0.
     651        END DO
     652      END IF
    629653      ok_sync=.TRUE.
    630654      IF (nqmax .LT. 2) THEN
    … …  
    644668c
    645669       IF (debut) THEN
    646 
     670C
     671         IF (if_ebil.ge.1) d_h_vcol_phy=0.
    647672c
    648673c appel a la lecture du run.def physique
    … …  
    15261551     .                "inst(X)", zsto,zout)
    15271552c
     1553         CALL histdef(nid_ins, "cdrm", "Momentum drag coef.", "-",
     1554     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     1555     .                "inst(X)", zsto,zout)
     1556c
     1557         CALL histdef(nid_ins, "cdrh", "Heat drag coef.", "-",
     1558     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     1559     .                "inst(X)", zsto,zout)
     1560c
    15281561         CALL histdef(nid_ins, "rain", "Precipitation", "mm/day",
    15291562     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
    … …  
    17311764         q_seri(i,k)  = qx(i,k,ivap)
    17321765         ql_seri(i,k) = qx(i,k,iliq)
     1766         qs_seri(i,k) = 0.
    17331767      ENDDO
    17341768      ENDDO
    … …  
    17481782      ENDDO
    17491783      ENDIF
    1750 c
     1784C
     1785      IF (if_ebil.ge.1) THEN
     1786        DO i = 1, klon
     1787          ztsol(i) = 0.
     1788        ENDDO
     1789        DO i = 1, klon
     1790          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
     1791        ENDDO
     1792        ztit='after dynamic'
     1793        CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
     1794     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     1795     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     1796C     Comme les tendances de la physique sont ajoute dans la dynamique,
     1797C     on devrait avoir que la variation d'entalpie par la dynamique
     1798C     est egale a la variation de la physique au pas de temps precedent.
     1799C     Donc la somme de ces 2 variations devrait etre nulle.
     1800        call diagphy(paire,ztit,ip_ebil
     1801     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     1802     e      , zero_v, zero_v, zero_v, ztsol
     1803     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
     1804     s      , fs_bound, fq_bound )
     1805      END IF
     1806
    17511807c Diagnostiquer la tendance dynamique
    17521808c
    … …  
    18121868      ENDDO
    18131869c
     1870      IF (if_ebil.ge.2) THEN
     1871        ztit='after reevap'
     1872        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
     1873     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     1874     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     1875         call diagphy(paire,ztit,ip_ebil
     1876     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     1877     e      , zero_v, zero_v, zero_v, ztsol
     1878     e      , d_h_vcol, d_qt, d_ec
     1879     s      , fs_bound, fq_bound )
     1880C
     1881      END IF
     1882C
     1883c
    18141884c Appeler la diffusion verticale (programme de couche limite)
    18151885c
    … …  
    18251895      DO nsrf = 1, nbsrf
    18261896      DO i = 1, klon
    1827          frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
     1897c$$$         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001
     1898        frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
    18281899      ENDDO
    18291900      ENDDO
    … …  
    18991970      ENDDO
    19001971      ENDDO
     1972c
     1973      IF (if_ebil.ge.2) THEN
     1974        ztit='after clmain'
     1975        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
     1976     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     1977     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     1978         call diagphy(paire,ztit,ip_ebil
     1979     e      , zero_v, zero_v, zero_v, zero_v, sens
     1980     e      , evap  , zero_v, zero_v, ztsol
     1981     e      , d_h_vcol, d_qt, d_ec
     1982     s      , fs_bound, fq_bound )
     1983      END IF
     1984C
    19011985c
    19021986c Incrementer la temperature du sol
    … …  
    20722156        ENDDO
    20732157      ENDDO
     2158c
     2159      IF (if_ebil.ge.2) THEN
     2160        ztit='after convect'
     2161        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
     2162     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     2163     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2164         call diagphy(paire,ztit,ip_ebil
     2165     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     2166     e      , zero_v, rain_con, snow_con, ztsol
     2167     e      , d_h_vcol, d_qt, d_ec
     2168     s      , fs_bound, fq_bound )
     2169      END IF
     2170C
    20742171      IF (check) THEN
    20752172          za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
    … …  
    21302227      ENDDO
    21312228      ENDDO
     2229c
     2230      IF (if_ebil.ge.2) THEN
     2231        ztit='after dry_adjust'
     2232        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
     2233     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     2234     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2235      END IF
    21322236
    21332237
    … …  
    21682272c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
    21692273c   relaxation des ratqs
    2170          facttemps=exp(-pdtphys/1.e4)
    2171          ratqs(:,:)=max(ratqs(:,:)*facttemps,ratqss(:,:))
     2274c         facttemps=exp(-pdtphys/1.e4)
     2275         facteur=exp(-pdtphys*facttemps)
     2276         ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
    21722277         ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
    21732278c         print*,'calcul des ratqs fini'
    … …  
    22142319      ENDIF
    22152320c
     2321      IF (if_ebil.ge.2) THEN
     2322        ztit='after fisrt'
     2323        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
     2324     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     2325     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2326        call diagphy(paire,ztit,ip_ebil
     2327     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     2328     e      , zero_v, rain_lsc, snow_lsc, ztsol
     2329     e      , d_h_vcol, d_qt, d_ec
     2330     s      , fs_bound, fq_bound )
     2331      END IF
     2332c
    22162333c-------------------------------------------------------------------
    22172334c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
    … …  
    22392356c  convection et du calcul du pas de temps précédent diminué d'un facteur
    22402357c  facttemps
    2241       facttemps=pdtphys/1.e4
     2358c      facttemps=pdtphys/1.e4
     2359      facteur = pdtphys *facttemps
    22422360      do k=1,klev
    22432361         do i=1,klon
    2244             rnebcon(i,k)=rnebcon(i,k)*facttemps
     2362            rnebcon(i,k)=rnebcon(i,k)*facteur
    22452363            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
    22462364     s      then
    … …  
    22782396         snow_fall(i) = snow_con(i) + snow_lsc(i)
    22792397      ENDDO
     2398c
     2399      IF (if_ebil.ge.2) THEN
     2400        ztit="after diagcld"
     2401        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
     2402     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     2403     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2404      END IF
    22802405c
    22812406c Calculer l'humidite relative pour diagnostique
    … …  
    23562481      ENDDO
    23572482c
     2483      IF (if_ebil.ge.2) THEN
     2484        ztit='after rad'
     2485        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
     2486     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     2487     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2488        call diagphy(paire,ztit,ip_ebil
     2489     e      , topsw, toplw, solsw, sollw, zero_v
     2490     e      , zero_v, zero_v, zero_v, ztsol
     2491     e      , d_h_vcol, d_qt, d_ec
     2492     s      , fs_bound, fq_bound )
     2493      END IF
     2494c
     2495c
    23582496c Calculer l'hydrologie de la surface
    23592497c
    … …  
    24572595c
    24582596      ENDIF ! fin de test sur ok_orolf
     2597c
     2598      IF (if_ebil.ge.2) THEN
     2599        ztit='after orography'
     2600        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
     2601     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     2602     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2603      END IF
     2604c
    24592605c
    24602606cAA
    … …  
    25202666c
    25212667c
    2522 
     2668      IF (if_ebil.ge.1) THEN
     2669        ztit='after physic'
     2670        CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
     2671     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     2672     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2673C     Comme les tendances de la physique sont ajoute dans la dynamique,
     2674C     on devrait avoir que la variation d'entalpie par la dynamique
     2675C     est egale a la variation de la physique au pas de temps precedent.
     2676C     Donc la somme de ces 2 variations devrait etre nulle.
     2677        call diagphy(paire,ztit,ip_ebil
     2678     e      , topsw, toplw, solsw, sollw, sens
     2679     e      , evap, rain_fall, snow_fall, ztsol
     2680     e      , d_h_vcol, d_qt, d_ec
     2681     s      , fs_bound, fq_bound )
     2682C
     2683      d_h_vcol_phy=d_h_vcol
     2684C
     2685      END IF
     2686C
    25232687      IF (ok_journe) THEN
    25242688c
    … …  
    32723436      CALL histwrite(nid_ins,"snow",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
    32733437
     3438      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragm,zx_tmp_2d)
     3439      CALL histwrite(nid_ins,"cdrm",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
     3440c
     3441      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragh,zx_tmp_2d)
     3442      CALL histwrite(nid_ins,"cdrh",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
    32743443c
    32753444      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
Note: See TracChangeset for help on using the changeset viewer.