Changeset 1411 for LMDZ4/trunk/libf


Ignore:
Timestamp:
Jul 9, 2010, 1:06:15 PM (14 years ago)
Author:
idelkadi
Message:

Nettoyage dans physiq.F des anciennes versions versions des options iflag_cldcon=5 ou 6.
Dans isrtilp.F, iflag_cldcon=5 : la version bi-gaussienne partout et iflag_cldcon=6 : bi-gaussiennes pour les thermiques et la lognormale ailleurs

Location:
LMDZ4/trunk/libf/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/phylmd/cloudth.F90

    r1403 r1411  
    142142! Calcul des écart-types pour s
    143143!-----------------------------------------------------------------------------------------------------------------
    144 !      alpha=0.5*(fraca(ind1,ind2)+fraca(ind1,ind2-1))
    145 !      sigma1s=(Max(1.2*fraca(ind1,ind2)*(sth-senv)**2,(senv/100)**2))**0.5
    146 !      sigma2s=((0.021/(fraca(ind1,ind2)+0.0)**0.5)*(sth-senv)**2+(sth/100)**2)**0.5
    147 !      sigma1s=(1.5**0.5)*(fraca(ind1,ind2)**0.65)*(sth-senv)+0.000045
    148 !      sigma2s=0.1265*(sth-senv)/(fraca(ind1,ind2)+0.0)**0.35+0.000045     
    149 
    150       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.00003
    151       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
    152 
    153 !      sigma1s=(1.5**0.5)*(alpha**0.65)*(sth-senv)+0.000045
    154 !      sigma2s=0.1265*(sth-senv)/alpha**0.35+0.000045
    155      
    156 !      sigma1s=(2.8**0.5)*(0.1**0.7)*(sth-senv)+0.00002
    157 !      sigma2s=((0.126/(0.1)**0.3)*(sth-senv))+0.00002
    158 
     144
     145      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
     146      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) 
    159147
    160148 
    … …  
    214202      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
    215203     
    216 !      if (zqenv(ind1).gt.0.005) then
    217       sigma1s=0.005*zqenv(ind1)
    218 !      else
    219 !      sigma1s=0.0005
    220 !      endif
    221 
    222 !      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
    223 
    224 
    225 
    226 !      sigma1s=0.00003
     204
     205      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
     206
    227207      sqrt2pi=sqrt(2.*pi)
    228208      xenv=senv/(sqrt(2.)*sigma1s)

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

    r1407 r1411  
    5050      REAL zthl(klon,klev)
    5151
     52      logical lognormale(klon)
     53
    5254cAA
    5355c Coeffients de fraction lessivee : pour OFF-LINE
    … …  
    140142      zdelq=0.0
    141143     
    142 !      print*,'CLOUDTH4 A. JAM'
     144      print*,'NUAGES4 A. JAM'
    143145      IF (appel1er) THEN
    144146c
    … …  
    289291     .                            /RG/dtime
    290292
    291 c pour la glace, on r�vapore toute la pr�ip dans la couche du dessous
     293c pour la glace, on r�vapore toute la pr�ip dans la couche du dessous
    292294c la glace venant de la couche du dessus est simplement dans la couche
    293295c du dessous.
    … …  
    343345c   zqn   : eau totale dans le nuage
    344346c   zcond : eau condensee moyenne dans la maille.
    345 c           on prend en compte le r�hauffement qui diminue la partie condensee
     347c           on prend en compte le r�hauffement qui diminue la partie condensee
    346348c
    347349c   Version avec les raqts
    … …  
    365367           enddo
    366368
    367               if (iflag_cldcon.eq.5) then
     369              if (iflag_cldcon>=5) then
    368370
    369371                 call cloudth(klon,klev,k,ztv,
    … …  
    377379                 enddo
    378380
    379               else
    380 
     381              endif
     382
     383! Pour iflag_cldcon<=4, on prend toujours la lognormale
     384! Dans le cas iflag_cldcon=5, on prend systématiquement la bi-gaussienne
     385! Dans le cas iflagÃ_cldcon=6, on prend la lognormale en absence des thermiques
     386
     387            lognormale(:)=
     388     .      iflag_cldcon<=4.or.(iflag_cldcon==6.and.fraca(:,k)<1.e-10)
    381389            do i=1,klon
     390            if (lognormale(i)) then
    382391            zpdf_sig(i)=ratqs(i,k)*zq(i)
    383392            zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
    … …  
    391400            zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
    392401            zpdf_e2(i)=1.-erf(zpdf_e2(i))
     402            endif
     403            enddo
     404
     405            do i=1,klon
     406            if (lognormale(i)) then
    393407            if (zpdf_e1(i).lt.1.e-10) then
    394408               rneb(i,k)=0.
    … …  
    398412               zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
    399413            endif
     414            endif
    400415           
    401416           enddo
    402417
    403          endif ! iflag_cldcon
    404418
    405419        endif ! iflag_pdf

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

    r1403 r1411  
    11! $Id$
     2!
    23c#define IO_DEBUG
    34
    … …  
    261262      CHARACTER*4 bb2
    262263      CHARACTER*2 bb3
     264c
    263265
    264266      real twriteSTD(klon,nlevSTD,nfiles)
    … …  
    863865      REAL rflag(klon)          ! flag fonctionnement de convect
    864866      INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
    865 
    866867c -- convect43:
    867868      INTEGER ntra              ! nb traceurs pour convect4.3
    … …  
    12321233      call phys_state_var_init(read_climoz)
    12331234      call phys_output_var_init
     1235
    12341236      print*, '================================================='
    12351237cIM for NMC files
    … …  
    12511253c         pmflxr=0.
    12521254c         pmflxs=0.
    1253 
    12541255        itau_con=0
    12551256        first=.false.
    … …  
    14231424           ema_pcb(i)  = 0.
    14241425           ema_pct(i)  = 0.
    1425 c          ema_workcbmf(i) = 0.
     1426           ema_workcbmf(i) = 0.
    14261427          ENDDO
    14271428cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
    … …  
    14481449c================================================================================
    14491450
    1450          ENDIF !debut
     1451         ENDIF
    14511452
    14521453           DO i=1,klon
    … …  
    21702171     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
    21712172     .        rain_con, snow_con, ibas_con, itop_con, sigd,
    2172      .        ema_cbmf,upwd,dnwd,dnwd0,
     2173     .        upwd,dnwd,dnwd0,
    21732174     .        Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
    21742175     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
    … …  
    21892190
    21902191          ELSE ! ok_cvl
    2191 
    21922192c MAF conema3 ne contient pas les traceurs
    21932193          CALL conema3 (dtime,
    … …  
    22292229
    22302230          DO i = 1, klon
    2231             ema_pcb(i)  = paprs(i,ibas_con(i))
     2231            ema_pcb(i)  = pbase(i)
    22322232          ENDDO
    22332233          DO i = 1, klon
     2234
    22342235! L'idicage de itop_con peut cacher un pb potentiel
    22352236! FH sous la dictee de JYG, CR
    … …  
    22422243              endif
    22432244            endif
    2244           ENDDO     
     2245          ENDDO
     2246          DO i = 1, klon
     2247            ema_cbmf(i) = ema_workcbmf(i)
     2248          ENDDO     
    22452249      ELSE IF (iflag_con.eq.0) THEN
    22462250          write(lunout,*) 'On n appelle pas la convection'
    … …  
    24592463c  ==============
    24602464
    2461 ! Dans le cas où on active les thermiques, on fait partir l'ajustement
     2465! Dans le cas où on active les thermiques, on fait partir l'ajustement
    24622466! a partir du sommet des thermiques.
    24632467! Dans le cas contraire, on demarre au niveau 1.
    … …  
    25442548         if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
    25452549
    2546 c-----------------------------------------------------------------------
    2547 c   par calcul direct de l'ecart-type
    2548 c-----------------------------------------------------------------------
    2549 
    2550       else if (iflag_cldcon>=5) then
    2551          wmax_th(:)=0.
    2552          zmax_th(:)=0.
    2553          do k=1,klev
    2554             do i=1,klon
    2555                wmax_th(i)=max(wmax_th(i),zw2(i,k))
    2556                if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
    2557             enddo
    2558          enddo
    2559          tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
    2560          if(prt_level.ge.9)
    2561      &     write(lunout,*)'TAU TH OK ',
    2562      &     tau_overturning_th(1),detr_therm(1,3)
    2563 
    2564 c On impose que l'air autour de la fraction couverte par le thermique
    2565 c plus son air detraine durant tau_overturning_th soit superieur
    2566 c a 0.1 q_seri
    2567          zz=0.1
    2568          do k=1,klev
    2569             do i=1,klon
    2570                lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
    2571      s         tau_overturning_th(i)*detr_therm(i,k)
    2572      s         *rg/(paprs(i,k)-paprs(i,k+1))
    2573                znum=(1.-zz)*q_seri(i,k)
    2574                zden=zqasc(i,k)-zz*q_seri(i,k)
    2575                if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
    2576                lambda_th(i,k)=min(lambda_th(i,k),0.9)
    2577             enddo
    2578          enddo
    2579 
    2580          if(iflag_cldcon==5) then
    2581             do k=1,klev
    2582                do i=1,klon
    2583                   ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
    2584      s            abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
    2585                enddo
    2586             enddo
    2587          else if(iflag_cldcon==6) then
    2588             do k=1,klev
    2589                do i=1,klon
    2590                   ratqsc(i,k)=sqrt(lambda_th(i,k))*
    2591      s            (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
    2592                enddo
    2593             enddo
    2594          endif
    2595 
    25962550      endif
    25972551
    … …  
    26572611
    26582612      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
    2659      s    .or.iflag_cldcon.ge.4) then
     2613     s    .or.iflag_cldcon.eq.4) then
    26602614
    26612615! On ajoute une constante au ratqsc*2 pour tenir compte de
    … …  
    26922646      endif
    26932647
     2648      print*,'PHSYIQ NUAGES4'
    26942649
    26952650c
    … …  
    28692824
    28702825c   On prend la somme des fractions nuageuses et des contenus en eau
    2871       cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
    2872       cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
     2826
     2827      if (iflag_cldcon>=5) then
     2828
     2829! Si on est sur un point touche par la convection profonde et pas
     2830! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
     2831! de la convection profonde.
     2832
     2833         print*,'TEST SCHEMA DE NUAGES '
     2834         do k=1,klev
     2835            do i=1,klon
     2836               if (ptconv(i,k).and. .not. ptconvth(i,k)) then
     2837                   cldfra(i,k)=rnebcon(i,k)
     2838                   cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
     2839               endif
     2840            enddo
     2841         enddo
     2842      else
     2843! Ancienne version
     2844         cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
     2845         cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
     2846      endif
    28732847
    28742848      ENDIF
    … …  
    33383312!     s        ref_liq,ref_ice
    33393313          call phys_cosp(itap,dtime,freq_cosp,
    3340      $                   ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP,
    3341      $                   ecrit_mth,ecrit_day,ecrit_hf,
    3342      $                   klon,klev,rlon,rlat,presnivs,overlap,
     3314     $                 ecrit_mth,ecrit_day,ecrit_hf,overlap,
     3315     $                   klon,klev,rlon,rlat,presnivs,
    33433316     $                   ref_liq,ref_ice,
    33443317     $                   pctsrf(:,is_ter)+pctsrf(:,is_lic),
    … …  
    33493322     $                   pmflxr(:,1:klev),pmflxs(:,1:klev),
    33503323     $                   mr_ozone,cldtau, cldemi)
    3351 
    33523324!     L          calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
    33533325!     L          cfaddbze,clcalipso2,dbze,cltlidarradar,
    … …  
    34973469      wwriteSTD(:,:,4)=wlevSTD(:,:)
    34983470c
    3499 cIM initialisation 5eme fichier de sortie
    35003471cIM ajoute 5eme niveau 170310 BEG
    35013472      twriteSTD(:,:,5)=tlevSTD(:,:)
    … …  
    36103581cIM global posePB#include "write_bilKP_ave.h"
    36113582c
    3612 
    36133583c Sauvegarder les valeurs de t et q a la fin de la physique:
    36143584c
    … …  
    36693639      DO k = 1, klev
    36703640      DO i = 1, klon
    3671 cJYG/IM theta en debut du pas de temps
    3672 cJYG/IM       theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
    3673 cJYG/IM theta en fin de pas de temps de physique
    3674         theta(i,k)=t_seri(i,k)*(100000./pplay(i,k))**(RD/RCPD)
     3641        theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
    36753642      ENDDO
    36763643      ENDDO
Note: See TracChangeset for help on using the changeset viewer.