Ignore:
Timestamp:
Jan 14, 2008, 1:03:39 PM (16 years ago)
Author:
Laurent Fairhead
Message:

Bascule de la physique du LMD vers la physique avec thermiques
LF

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/phylmd/physiq.F

    r864 r878  
    177177      REAL fm_therm(klon,klev+1)
    178178      REAL entr_therm(klon,klev)
     179      real,allocatable,save :: clwcon0th(:,:),rnebcon0th(:,:)
     180c$OMP THREADPRIVATE(clwcon0th,rnebcon0th)
     181
     182      real weak_inversion(klon),dthmin(klon)
     183      real seuil_inversion
     184      save seuil_inversion
     185c$OMP THREADPRIVATE(seuil_inversion)
     186      integer iflag_ratqs
     187      save iflag_ratqs
     188c$OMP THREADPRIVATE(iflag_ratqs)
     189
     190      integer lmax_th(klon)
     191      integer limbas(klon)
     192      real ratqscth(klon,klev)
     193      real ratqsdiff(klon,klev)
     194      real zqsatth(klon,klev)
     195
    179196c======================================================================
    180197c
    … …  
    929946      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
    930947      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
     948
     949      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: pbl_tke ! turb kinetic energy
     950c   !$OMP THREADPRIVATE(pbl_tke)
     951
    931952c
    932953      REAL zxfluxt(klon, klev)
    … …  
    10501071c$OMP THREADPRIVATE(d_u_con,d_v_con)
    10511072      REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
     1073      REAL d_t_ajsb(klon,klev), d_q_ajsb(klon,klev)
    10521074      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
    10531075      REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
    … …  
    14711493#endif
    14721494      allocate( clwcon0(klon,klev),rnebcon0(klon,klev))
     1495      allocate( clwcon0th(klon,klev),rnebcon0th(klon,klev))
    14731496      allocate( tau_ae(klon,klev,2), piz_ae(klon,klev,2))
    14741497      allocate( cg_ae(klon,klev,2))
    … …  
    15071530         CALL suphec ! initialiser constantes et parametres phys.
    15081531      ENDIF
     1532
     1533      print*,'CONVERGENCE PHYSIQUE THERM 1 '
    15091534
    15101535
    … …  
    15671592c
    15681593         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
    1569      .                  ok_instan, ok_hf,
     1594     .                  ok_instan, ok_hf, seuil_inversion,
    15701595     .                  fact_cldcon, facttemps,ok_newmicro,
    1571 cIM  .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
    1572      .                  iflag_cldcon,ratqsbas,ratqshaut,
     1596     .                  iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
    15731597     .                  ok_ade, ok_aie,
    15741598     .                  bl95_b0, bl95_b1,
    … …  
    15811605         itap    = 0
    15821606         itaprad = 0
     1607
     1608!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1609!! Un petit travail à faire ici.
     1610!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1611
     1612         if (iflag_pbl>1) then
     1613             PRINT*, "Using method MELLOR&YAMADA"
     1614         endif
     1615          ! NB! pbl_tke could/should be read and written from (re)startphy.nc
     1616          ALLOCATE(pbl_tke(klon,klev+1,nbsrf)) 
     1617          pbl_tke(:,:,:) = 1.e-8
     1618
    15831619
    15841620         CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
    … …  
    15891625     .       radsol,clesphy0,
    15901626     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
    1591      .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon)
     1627     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
     1628     .       pbl_tke)
     1629
     1630
     1631
     1632
     1633!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1634
    15921635
    15931636       DO i=1,klon
    … …  
    16201663     .                        pdtphys
    16211664            abort_message='Pas physique n est pas correct '
    1622             call abort_gcm(modname,abort_message,1)
     1665!           call abort_gcm(modname,abort_message,1)
     1666            dtime=pdtphys
    16231667         ENDIF
    16241668         IF (nlon .NE. klon) THEN
    … …  
    16701714         ENDIF
    16711715
     1716           rugoro=0.
    16721717c34EK
    16731718         IF (ok_orodr) THEN
    1674            DO i=1,klon
    1675              rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
    1676            ENDDO
     1719
     1720           rugoro=0.
     1721
     1722!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1723! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
     1724! justement quand ok_orodr = false.
     1725! ce rugoro est utilise par la couche limite et fait double emploi
     1726! avec les paramétrisations spécifiques de Francois Lott.
     1727!           DO i=1,klon
     1728!             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
     1729!           ENDDO
     1730!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1731
    16771732           CALL SUGWD(klon,klev,paprs,pplay)
    16781733           DO i=1,klon
    … …  
    20882143     d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
    20892144     -     dsens,     devap,     zxsnow,
    2090      -     zxfluxt,   zxfluxq,   q2m,     fluxq )
     2145     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
    20912146c
    20922147c
    … …  
    22602315          DO i = 1, klon
    22612316            ema_pct(i)  = paprs(i,itop_con(i))
     2317            if (itop_con(i).gt.klev-3) then
     2318               print*,'La convection monte trop haut '
     2319               print*,'itop_con(,',i,',)=',itop_con(i)
     2320            endif
    22622321          ENDDO
    22632322          DO i = 1, klon
    … …  
    23432402c===================================================================
    23442403c
     2404       call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
     2405     s ,seuil_inversion,weak_inversion,dthmin)
     2406
     2407
     2408
     2409      d_t_ajsb(:,:)=0.
     2410      d_q_ajsb(:,:)=0.
    23452411      d_t_ajs(:,:)=0.
    23462412      d_u_ajs(:,:)=0.
    23472413      d_v_ajs(:,:)=0.
    23482414      d_q_ajs(:,:)=0.
     2415      clwcon0th(:,:)=0.
    23492416      fm_therm(:,:)=0.
    23502417      entr_therm(:,:)=0.
    … …  
    23572424c  ====
    23582425         IF(prt_level>9)WRITE(lunout,*)'pas de convection'
    2359       else if(iflag_thermals.eq.0) then
    2360 
    2361 c  Ajustement sec
    2362 c  ==============
    2363          IF(prt_level>9)WRITE(lunout,*)'ajsec'
    2364          CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
    2365          t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
    2366          q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
     2426
     2427
    23672428      else
     2429
    23682430c  Thermiques
    23692431c  ==========
    23702432         IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
    23712433     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
     2434         print*,'JUSTE AVANT , iflag_thermals='
     2435     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
     2436
     2437
     2438         if (iflag_thermals.gt.1) then
    23722439         call calltherm(pdtphys
    2373      s      ,pplay,paprs,pphi
    2374      s      ,u_seri,v_seri,t_seri,q_seri
     2440     s      ,pplay,paprs,pphi,weak_inversion
     2441     s      ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
    23752442     s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
    2376      s      ,fm_therm,entr_therm)
     2443     s      ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth,
     2444     s       ratqsdiff,zqsatth)
     2445         endif
     2446
     2447!        call calltherm(pdtphys
     2448!    s      ,pplay,paprs,pphi
     2449!    s      ,u_seri,v_seri,t_seri,q_seri
     2450!    s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
     2451!    s      ,fm_therm,entr_therm)
     2452
     2453c  Ajustement sec
     2454c  ==============
     2455
     2456! Dans le cas où on active les thermiques, on fait partir l'ajustement
     2457! a partir du sommet des thermiques.
     2458! Dans le cas contraire, on demarre au niveau 1.
     2459
     2460         if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
     2461
     2462         if(iflag_thermals.eq.0) then
     2463            IF(prt_level>9)WRITE(lunout,*)'ajsec'
     2464            limbas(:)=1
     2465         else
     2466            limbas(:)=lmax_th(:)
     2467         endif
     2468 
     2469! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
     2470! pour des test de convergence numerique.
     2471! Le nouveau ajsec est a priori mieux, meme pour le cas
     2472! iflag_thermals = 0 (l'ancienne version peut faire des tendances
     2473! non nulles numeriquement pour des mailles non concernees.
     2474
     2475         if (iflag_thermals.eq.0) then
     2476            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
     2477     s      , d_t_ajsb, d_q_ajsb)
     2478         else
     2479            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
     2480     s      , d_t_ajsb, d_q_ajsb)
     2481         endif
     2482
     2483         t_seri(:,:) = t_seri(:,:) + d_t_ajsb(:,:)
     2484         q_seri(:,:) = q_seri(:,:) + d_q_ajsb(:,:)
     2485         d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
     2486         d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
     2487
     2488         endif
     2489
    23772490      endif
    23782491c
    … …  
    24062519         enddo
    24072520         enddo
     2521      else if (iflag_cldcon.eq.4) then
     2522         ptconvth(:,:)=.false.
     2523         ratqsc(:,:)=0.
     2524         print*,'avant clouds_gno thermique'
     2525         call clouds_gno
     2526     s   (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
     2527         print*,' CLOUDS_GNO OK'
    24082528      endif
    24092529
    24102530c   ratqs stables
    24112531c   -------------
    2412       do k=1,klev
    2413 cIM RAJOUT boucle do=i
    2414        do i=1, klon
    2415 cIM         ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
    2416 cIM     s   min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
    2417          ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
    2418      s   min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
    2419 cIM      print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
    2420        enddo
    2421       enddo
     2532
     2533      if (iflag_ratqs.eq.0) then
     2534
     2535! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
     2536         do k=1,klev
     2537            do i=1, klon
     2538               ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
     2539     s         min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
     2540            enddo
     2541         enddo
     2542
     2543! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
     2544! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
     2545! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
     2546! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
     2547! Il s'agit de differents tests dans la phase de reglage du modele
     2548! avec thermiques.
     2549
     2550      else if (iflag_ratqs.eq.1) then
     2551
     2552         do k=1,klev
     2553            do i=1, klon
     2554               if (pplay(i,k).ge.60000.) then
     2555                  ratqss(i,k)=ratqsbas
     2556               else if ((pplay(i,k).ge.30000.).and.
     2557     s            (pplay(i,k).lt.60000.)) then
     2558                  ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
     2559     s            (60000.-pplay(i,k))/(60000.-30000.)
     2560               else
     2561                  ratqss(i,k)=ratqshaut
     2562               endif
     2563            enddo
     2564         enddo
     2565
     2566      else
     2567
     2568         do k=1,klev
     2569            do i=1, klon
     2570               if (pplay(i,k).ge.60000.) then
     2571                  ratqss(i,k)=ratqsbas
     2572     s            *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
     2573               else if ((pplay(i,k).ge.30000.).and.
     2574     s             (pplay(i,k).lt.60000.)) then
     2575                    ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
     2576     s              (60000.-pplay(i,k))/(60000.-30000.)
     2577               else
     2578                    ratqss(i,k)=ratqshaut
     2579               endif
     2580            enddo
     2581         enddo
     2582      endif
     2583
     2584
    24222585
    24232586
    24242587c  ratqs final
    24252588c  -----------
    2426       if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
    2427 c   les ratqs sont une conbinaison de ratqss et ratqsc
    2428 c   ratqs final
    2429 c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
    2430 c   relaxation des ratqs
    2431 c         facttemps=exp(-pdtphys/1.e4)
    2432          facteur=exp(-pdtphys*facttemps)
    2433          ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
    2434          ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
    2435 c         print*,'calcul des ratqs fini'
     2589
     2590      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
     2591     s    .or.iflag_cldcon.eq.4) then
     2592
     2593! On ajoute une constante au ratqsc*2 pour tenir compte de
     2594! fluctuations turbulentes de petite echelle
     2595
     2596         do k=1,klev
     2597            do i=1,klon
     2598               if ((fm_therm(i,k).gt.1.e-10)) then
     2599                  ratqsc(i,k)=sqrt(ratqsc(i,k)**2+0.05**2)
     2600               endif
     2601            enddo
     2602         enddo
     2603
     2604!   les ratqs sont une conbinaison de ratqss et ratqsc
     2605!   1800s, en dur pour le moment, est le temps de
     2606!   relaxation des ratqs
     2607
     2608         facteur=exp(-pdtphys/1800.)
     2609
     2610         print*,'WARNING ratqs a revoir '
     2611         ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
     2612         ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
     2613
    24362614      else
    2437 c   on ne prend que le ratqs stable pour fisrtilp
     2615!   on ne prend que le ratqs stable pour fisrtilp
    24382616         ratqs(:,:)=ratqss(:,:)
    24392617      endif
    … …  
    25052683cIM cf FH
    25062684c     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
    2507        IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
     2685      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
    25082686       snow_tiedtke=0.
    25092687c     print*,'avant calcul de la pseudo precip '
    … …  
    25422720      ENDDO
    25432721
    2544       ELSE IF (iflag_cldcon.eq.3) THEN
     2722      ELSE IF (iflag_cldcon.ge.3) THEN
    25452723c  On prend pour les nuages convectifs le max du calcul de la
    25462724c  convection et du calcul du pas de temps precedent diminue d'un facteur
    25472725c  facttemps
    2548 c      facttemps=pdtphys/1.e4
    25492726      facteur = pdtphys *facttemps
    25502727      do k=1,klev
    … …  
    33143491     .      radsol,
    33153492     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
    3316      .      t_ancien, q_ancien, rnebcon, ratqs, clwcon)
     3493     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon,
     3494     .      pbl_tke)
    33173495      ENDIF
    33183496     
Note: See TracChangeset for help on using the changeset viewer.