Ignore:
Timestamp:
Jul 12, 2002, 12:27:22 PM (22 years ago)
Author:
lmdzadmin
Message:

Inclusion du nouveau schema de nuages de SB. FH
IM/LF

Location:
LMDZ.3.3/branches/rel-LF/libf/phylmd
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r230 r373  
    11c $Header$
    2       SUBROUTINE conema (dtime,paprs,pplay,t,q,u,v,tra,ntra,
     2      SUBROUTINE conema3 (dtime,paprs,pplay,t,q,u,v,tra,ntra,
    33     .             work1,work2,d_t,d_q,d_u,d_v,d_tra,
    44     .             rain, snow, kbas, ktop,
    55     .             upwd,dnwd,dnwdbis,bas,top,Ma,cape,tvp,rflag,
    6      .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
     6     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
     7     .             qcond_incld)
    78
    89      IMPLICIT none
    … …  
    5253#include "dimensions.h"
    5354#include "dimphy.h"
     55#include "conema3.h"
    5456      INTEGER i, l,m,itra
    5557      INTEGER ntra,ntrac !number of tracers; if no tracer transport
    … …  
    5860      REAL dtime
    5961c
    60 
     62      REAL d_t2(klon,klev), d_q2(klon,klev) ! sbl
     63      REAL d_u2(klon,klev), d_v2(klon,klev) ! sbl
     64      REAL em_d_t2(klev), em_d_q2(klev)     ! sbl   
     65      REAL em_d_u2(klev), em_d_v2(klev)     ! sbl   
    6166c
    6267      REAL paprs(klon,klev+1), pplay(klon,klev)
    … …  
    7378      REAL dplcldt(klon), dplcldr(klon)
    7479      INTEGER kbas(klon), ktop(klon)
     80
     81      REAL wd(klon)
     82      REAL qcond_incld(klon,klev)
    7583c
    7684      REAL em_t(klev)
    … …  
    9199      INTEGER em_bas, em_top
    92100      SAVE em_bas, em_top
     101
     102      REAL em_wd
     103      REAL em_qcond(klev)
     104      REAL em_qcondc(klev)
    93105c
    94106      REAL zx_t, zx_qs, zdelta, zcor
    … …  
    115127      real em_pbase, em_bbase
    116128      integer iw,j,k,ix,iy
     129
     130c -- sb: pour schema nuages:
     131
     132       integer iflagcon
     133       integer em_ifc(klev)
     134     
     135       real em_pradj
     136       real em_cldf(klev), em_cldq(klev)
     137       real em_ftadj(klev), em_fradj(klev)
     138
     139       integer ifc(klon,klev)
     140       real pradj(klon)
     141       real cldf(klon,klev), cldq(klon,klev)
     142       real ftadj(klon,klev), fqadj(klon,klev)
     143
     144c sb --
     145
    117146ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    118147c
    … …  
    121150#include "FCTTRE.h"
    122151 
     152      qcond_incld(:,:) = 0.
    123153c
    124154c$$$      print*,'debut conema'
    … …  
    190220ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    191221c
    192       CALL convect(dtime, em_t, em_q, em_qs,em_u ,em_v ,
     222
     223c     print*,'avant convect i=',i
     224      CALL convect3(dtime,epmax,ok_adj_ema,
     225     .              em_t, em_q, em_qs,em_u ,em_v ,
    193226     .              em_tra, em_p, em_ph,
    194227     .              klev, klev+1, klev-1,ntra, dtime, iflag,
    … …  
    197230     .              em_bas, em_top,em_upwd, em_dnwd, em_dnwdbis,
    198231     .              em_work1, em_work2,emmip,emMke,emMa,Ment,
    199 c SB 11sept98     .              Qent,TPS,TLS,SIJ)
    200 c 19oct98     .              Qent,TPS,TLS,SIJ,em_CAPE,em_TVP)
    201232     .  Qent,TPS,TLS,SIJ,em_CAPE,em_TVP,em_pbase,em_bbase,
    202      .  em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr)
     233     .  em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr, ! sbl
     234     .  em_d_t2,em_d_q2,em_d_u2,em_d_v2,em_wd,em_qcond,em_qcondc)!sbl
     235c     print*,'apres convect '
    203236c
    204237ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     238c
     239c -- sb: Appel schema statistique de nuages couple a la convection
     240c (Bony et Emanuel 2001):
     241
     242c -- creer cvthermo.h qui contiendra les cstes thermo de LMDZ:
     243
     244        iflagcon = 3
     245c       CALL cv_thermo(iflagcon)
     246
     247c -- appel schema de nuages:
     248
     249c       CALL CLOUDS_SUB_LS(klev,em_q,em_qs,em_t
     250c    i          ,em_p,em_ph,dtime,em_qcondc
     251c    o          ,em_cldf,em_cldq,em_pradj,em_ftadj,em_fradj,em_ifc)
     252
     253        do k = 1, klev
     254         cldf(i,k)  = em_cldf(k)  ! cloud fraction (0-1)
     255         cldq(i,k)  = em_cldq(k)  ! in-cloud water content (kg/kg)
     256         ftadj(i,k) = em_ftadj(k) ! (dT/dt)_{LS adj} (K/s)
     257         fqadj(i,k) = em_fradj(k) ! (dq/dt)_{LS adj} (kg/kg/s)
     258         ifc(i,k)   = em_ifc(k)   ! flag convergence clouds_gno (1 ou 2)
     259        enddo
     260        pradj(i) = em_pradj       ! precip from LS supersat adj (mm/day)
     261
     262c sb --
    205263c
    206264c SB:
    … …  
    251309      snow(i) = 0.0
    252310      cape(i) = em_CAPE
     311      wd(i) = em_wd
    253312      rflag(i) = float(iflag)
    254313c SB      kbas(i) = em_bas
    … …  
    257316      dplcldr(i) = em_dplcldr
    258317      DO l = 1, klev
     318         d_t2(i,l) = dtime * em_d_t2(l)
     319         d_q2(i,l) = dtime * em_d_q2(l)
     320         d_u2(i,l) = dtime * em_d_u2(l)
     321         d_v2(i,l) = dtime * em_d_v2(l)
     322
    259323         d_t(i,l) = dtime * em_d_t(l)
    260324         d_q(i,l) = dtime * em_d_q(l)
    … …  
    273337         dtvpdt1(i,l) = em_dtvpdt1(l)
    274338         dtvpdq1(i,l) = em_dtvpdq1(l)
     339
     340         if (iflag_clw.eq.0) then
     341            qcond_incld(i,l) = em_qcondc(l)
     342         else if (iflag_clw.eq.1) then
     343            qcond_incld(i,l) = em_qcond(l)
     344         endif
    275345      ENDDO
    276346  999 CONTINUE
     347
     348c   On calcule une eau liquide diagnostique en fonction de la
     349c  precip.
     350      if ( iflag_clw.eq.2 ) then
     351      do l=1,klev
     352         do i=1,klon
     353            if (ktop(i)-kbas(i).gt.0.and.
     354     s         l.ge.kbas(i).and.l.le.ktop(i)) then
     355               qcond_incld(i,l)=rain(i)*8.e4
     356c    s         *(pplay(i,l      )-paprs(i,ktop(i)+1))
     357     s         /(pplay(i,kbas(i))-pplay(i,ktop(i)))
     358c    s         **2
     359            else
     360               qcond_incld(i,l)=0.
     361            endif
     362         enddo
     363         print*,'l=',l,',   qcond_incld=',qcond_incld(1,l)
     364      enddo
     365      endif
    277366 
    278367

  • LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90

    r241 r373  
    33!
    44
    5   subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan)
     5  subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &
     6 &                     fact_cldcon, facttemps,ok_newmicro,iflag_cldcon, &
     7 &                     ratqsbas,ratqshaut)
    68
    79   use IOIPSL
    810   implicit none
    911
     12#include "conema3.h"
     13#include "fisrtilp.inc"
    1014!
    1115! Configuration de la "physique" de LMDZ a l'aide de la fonction
    … …  
    2630! Sortie:
    2731  character (len = 6)  :: ocean
    28   logical              :: ok_veget
     32  logical              :: ok_veget, ok_newmicro
    2933  logical              :: ok_journe, ok_mensuel, ok_instan
     34  real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut
     35  integer              :: iflag_cldcon
    3036! Local
    3137  integer              :: numout = 6
    … …  
    7884  ok_instan = .false.
    7985  call getin('OK_instan', ok_instan)
     86!!
     87!! KE
     88!
     89!Config Key  = epmax
     90!Config Desc = Efficacite precip
     91!Config Def  = 0.993
     92!Config Help =
     93!
     94  epmax = .993
     95  call getin('epmax', epmax)
     96!
     97!Config Key  = ok_adj_ema
     98!Config Desc = 
     99!Config Def  = false
     100!Config Help =
     101!
     102  ok_adj_ema = .false.
     103  call getin('ok_adj_ema',ok_adj_ema)
     104!
     105!Config Key  = iflag_clw
     106!Config Desc = 
     107!Config Def  = 0
     108!Config Help =
     109!
     110  iflag_clw = 0
     111  call getin('iflag_clw',iflag_clw)
     112!
     113!Config Key  = cld_lc_lsc
     114!Config Desc = 
     115!Config Def  = 2.6e-4
     116!Config Help =
     117!
     118  cld_lc_lsc = 2.6e-4
     119  call getin('cld_lc_lsc',cld_lc_lsc)
     120!
     121!Config Key  = cld_lc_con
     122!Config Desc = 
     123!Config Def  = 2.6e-4
     124!Config Help =
     125!
     126  cld_lc_con = 2.6e-4
     127  call getin('cld_lc_con',cld_lc_con)
     128!
     129!Config Key  = cld_tau_lsc
     130!Config Desc = 
     131!Config Def  = 3600.
     132!Config Help =
     133!
     134  cld_tau_lsc = 3600.
     135  call getin('cld_tau_lsc',cld_tau_lsc)
     136!
     137!Config Key  = cld_tau_con
     138!Config Desc = 
     139!Config Def  = 3600.
     140!Config Help =
     141!
     142  cld_tau_con = 3600.
     143  call getin('cld_tau_con',cld_tau_con)
     144!
     145!Config Key  = ffallv_lsc
     146!Config Desc = 
     147!Config Def  = 1.
     148!Config Help =
     149!
     150  ffallv_lsc = 1.
     151  call getin('ffallv_lsc',ffallv_lsc)
     152!
     153!Config Key  = ffallv_con
     154!Config Desc = 
     155!Config Def  = 1.
     156!Config Help =
     157!
     158  ffallv_con = 1.
     159  call getin('ffallv_con',ffallv_con)
     160!
     161!Config Key  = coef_eva
     162!Config Desc = 
     163!Config Def  = 2.e-5
     164!Config Help =
     165!
     166  coef_eva = 2.e-5
     167  call getin('coef_eva',coef_eva)
     168!
     169!Config Key  = reevap_ice
     170!Config Desc = 
     171!Config Def  = .false.
     172!Config Help =
     173!
     174  reevap_ice = .false.
     175  call getin('reevap_ice',reevap_ice)
     176!
     177!Config Key  = iflag_cldcon
     178!Config Desc = 
     179!Config Def  = 1
     180!Config Help =
     181!
     182  iflag_cldcon = 1
     183  call getin('iflag_cldcon',iflag_cldcon)
     184
     185!
     186!Config Key  = iflag_pdf
     187!Config Desc = 
     188!Config Def  = 0
     189!Config Help =
     190!
     191  iflag_pdf = 0
     192  call getin('iflag_pdf',iflag_pdf)
     193!
     194!Config Key  = fact_cldcon
     195!Config Desc = 
     196!Config Def  = 0.375
     197!Config Help =
     198!
     199  fact_cldcon = 0.375
     200  call getin('fact_cldcon',fact_cldcon)
     201
     202!
     203!Config Key  = facttemps
     204!Config Desc = 
     205!Config Def  = 1.e-4
     206!Config Help =
     207!
     208  facttemps = 1.e-4
     209  call getin('facttemps',facttemps)
     210
     211!
     212!Config Key  = ok_newmicro
     213!Config Desc = 
     214!Config Def  = .true.
     215!Config Help =
     216!
     217  ok_newmicro = .true.
     218  call getin('ok_newmicro',ok_newmicro)
     219!
     220!Config Key  = ratqsbas
     221!Config Desc = 
     222!Config Def  = 0.01
     223!Config Help =
     224!
     225  ratqsbas = 0.01
     226  call getin('ratqsbas',ratqsbas)
     227!
     228!Config Key  = ratqshaut
     229!Config Desc = 
     230!Config Def  = 0.3
     231!Config Help =
     232!
     233  ratqshaut = 0.3
     234  call getin('ratqshaut',ratqshaut)
    80235
    81236!
    … …  
    90245  write(numout,*)' Sortie mensuelle = ', ok_mensuel
    91246  write(numout,*)' Sortie instantanee = ', ok_instan
    92 
     247  write(numout,*)' epmax = ', epmax
     248  write(numout,*)' ok_adj_ema = ', ok_adj_ema
     249  write(numout,*)' iflag_clw = ', iflag_clw
     250  write(numout,*)' cld_lc_lsc = ', cld_lc_lsc
     251  write(numout,*)' cld_lc_con = ', cld_lc_con
     252  write(numout,*)' cld_tau_lsc = ', cld_tau_lsc
     253  write(numout,*)' cld_tau_con = ', cld_tau_con
     254  write(numout,*)' ffallv_lsc = ', ffallv_lsc
     255  write(numout,*)' ffallv_con = ', ffallv_con
     256  write(numout,*)' coef_eva = ', coef_eva
     257  write(numout,*)' reevap_ice = ', reevap_ice
     258  write(numout,*)' iflag_pdf = ', iflag_pdf
     259  write(numout,*)' iflag_cldcon = ', iflag_cldcon
     260  write(numout,*)' fact_cldcon = ', fact_cldcon
     261  write(numout,*)' facttemps = ', facttemps
     262  write(numout,*)' ok_newmicro = ',ok_newmicro
     263  write(numout,*)' ratqsbas = ',ratqsbas
     264  write(numout,*)' ratqshaut = ',ratqshaut
    93265
    94266  return

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

    r230 r373  
    11c $Header$
    2       SUBROUTINE CONVECT
    3      *    (DTIME,  T1,   R1,   RS,    U,  V,  TRA,   P,     PH,
     2      SUBROUTINE CONVECT3       
     3     *    (DTIME,EPMAX,ok_adj,
     4     *     T1,   R1,   RS,    U,  V,  TRA,   P,     PH,
    45     *     ND,       NDP1,     NL, NTRA,  DELT,  IFLAG,
    56     *     FT, FR, FU,  FV,  FTRA,  PRECIP,
    67     *     icb,inb,   upwd,dnwd,dnwd0,SIG, W0,Mike,Mke,
    78     *     Ma,MENTS,QENTS,TPS,TLS,SIGIJ,CAPE,TVP,PBASE,BUOYBASE,
    8      *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR)
     9cccc     *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR)
     10     *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR,   ! sbl
     11     *     FT2,FR2,FU2,FV2,WD,QCOND,QCONDC)   ! sbl
    912C
    1013C    ***  THE PARAMETER NA SHOULD IN GENERAL EQUAL ND   ***
    … …  
    1922      integer NTRAC
    2023      PARAMETER (NTRAC=nqmx-2)
     24      REAL DELTAC              ! cld
     25      PARAMETER (DELTAC=0.01)  ! cld
    2126
    2227      INTEGER NENT(NA)
    … …  
    3439      REAL T(NA),RR(NA)
    3540C
     41      REAL FT2(ND),FR2(ND),FU2(ND),FV2(ND) ! added sbl
     42      REAL U1(ND),V1(ND) ! added sbl
     43C
    3644      REAL BUOY(NA)     !  Lifted parcel buoyancy
    3745      REAL DTVPDT1(ND),DTVPDQ1(ND)   ! Derivatives of parcel virtual
    … …  
    3947      REAL CLW_NEW(NA),QI(NA)
    4048C
    41       LOGICAL ICE_CONV
     49      REAL WD, BETAD ! for gust factor (sb)
     50      REAL QCONDC(ND)  ! interface cld param (sb)
     51      REAL QCOND(ND),NQCOND(NA),WA(NA),MAA(NA),SIGA(NA),AXC(NA) ! cld
     52C
     53      LOGICAL ICE_CONV,ok_adj
    4254      PARAMETER (ICE_CONV=.TRUE.)
    4355 
    … …  
    110122         FU(I)=0.0
    111123         FV(I)=0.0
     124
     125         FT2(I)=0.0
     126         FR2(I)=0.0
     127         FU2(I)=0.0
     128         FV2(I)=0.0
     129
    112130         DO 4 J=1,NTRA
    113131          FTRA(I,J)=0.0
    114132    4    CONTINUE
     133
     134         QCONDC(I)=0.0  ! cld
     135         QCOND(I)=0.0   ! cld
     136         NQCOND(I)=0.0  ! cld
     137
    115138         T(I)=T1(I)
    116139         RR(I)=R1(I)
     140         U1(I)=U(I) ! added sbl
     141         V1(I)=V(I) ! added sbl
    117142    5 CONTINUE
    118143      DO 7 I=1,NL
    … …  
    137162 
    138163      PRECIP=0.0
     164      WD=0.0 ! sb
    139165      IFLAG=1
    140166C
    … …  
    160186      SPFAC=0.15
    161187c sb:
    162       EPMAX=0.993 ! precip efficiency less than unity
     188c     EPMAX=0.993 ! precip efficiency less than unity
    163189c      EPMAX=1. ! precip efficiency less than unity
    164190C
    … …  
    196222C
    197223      NOPT=0
     224c!      NOPT=1 ! sbl
    198225C
    199226C     ***            PERFORM DRY ADIABATIC ADJUSTMENT            ***
    … …  
    202229C     ***                BOUNDARY LAYER SCHEME !!!               ***
    203230C
    204 c test sb      DO 30 I=NL-1,1,-1
    205 c test sb         JN=0
    206 c test sb         DO 10 J=I+1,NL
    207 c test sb   10    IF(TH(J).LT.TH(I))JN=J
    208 c test sb         IF(JN.EQ.0)GOTO 30
    209 c test sb         AHM=0.0
    210 c test sb         RM=0.0
    211 c test sb         UM=0.0
    212 c test sb         VM=0.0
    213 c test sb         DO K=1,NTRA
    214 c test sb          TRATM(K)=0.0
    215 c test sb         END DO
    216 c test sb         DO 15 J=I,JN
    217 c test sb          AHM=AHM+(CPD*(1.-RR(J))+RR(J)*CPV)*T(J)*(PH(J)-PH(J+1))
    218 c test sb          RM=RM+RR(J)*(PH(J)-PH(J+1))
    219 c test sb          UM=UM+U(J)*(PH(J)-PH(J+1))
    220 c test sb          VM=VM+V(J)*(PH(J)-PH(J+1))
    221 c test sb          DO K=1,NTRA
    222 c test sb           TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1))
    223 c test sb          END DO
    224 c test sb   15    CONTINUE
    225 c test sb         DPHINV=1./(PH(I)-PH(JN+1))
    226 c test sb         RM=RM*DPHINV
    227 c test sb         UM=UM*DPHINV
    228 c test sb         VM=VM*DPHINV
    229 c test sb         DO K=1,NTRA
    230 c test sb          TRATM(K)=TRATM(K)*DPHINV
    231 c test sb         END DO
    232 c test sb         A2=0.0
    233 c test sb         DO 20 J=I,JN
    234 c test sb            RR(J)=RM
    235 c test sb          U(J)=UM
    236 c test sb          V(J)=VM
    237 c test sb          DO K=1,NTRA
    238 c test sb           TRA(J,K)=TRATM(K)
    239 c test sb          END DO
    240 c test sb            RDCP=(RD*(1.-RR(J))+RR(J)*RV)/ (CPD*(1.-RR(J))+RR(J)*CPV)
    241 c test sb            X=(0.001*P(J))**RDCP
    242 c test sb            T(J)=X
    243 c test sb            A2=A2+(CPD*(1.-RR(J))+RR(J)*CPV)*X*(PH(J)-PH(J+1))
    244 c test sb   20    CONTINUE
    245 c test sb         DO 25 J=I,JN
    246 c test sb            TH(J)=AHM/A2
    247 c test sb            T(J)=T(J)*TH(J)
    248 c test sb   25    CONTINUE
    249 c test sb   30 CONTINUE
    250 C
    251 C   ***   RESET INPUT ARRAYS IF NOPT .NE. 0   ***
    252 C
    253       IF (NOPT.NE.0)THEN
     231      IF (ok_adj) THEN ! added sbl
     232
     233      DO 30 I=NL-1,1,-1
     234         JN=0
     235         DO 10 J=I+1,NL
     236   10    IF(TH(J).LT.TH(I))JN=J
     237         IF(JN.EQ.0)GOTO 30
     238         AHM=0.0
     239         RM=0.0
     240         UM=0.0
     241         VM=0.0
     242         DO K=1,NTRA
     243          TRATM(K)=0.0
     244         END DO
     245         DO 15 J=I,JN
     246          AHM=AHM+(CPD*(1.-RR(J))+RR(J)*CPV)*T(J)*(PH(J)-PH(J+1))
     247          RM=RM+RR(J)*(PH(J)-PH(J+1))
     248          UM=UM+U(J)*(PH(J)-PH(J+1))
     249          VM=VM+V(J)*(PH(J)-PH(J+1))
     250          DO K=1,NTRA
     251           TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1))
     252          END DO
     253   15    CONTINUE
     254         DPHINV=1./(PH(I)-PH(JN+1))
     255         RM=RM*DPHINV
     256         UM=UM*DPHINV
     257         VM=VM*DPHINV
     258         DO K=1,NTRA
     259          TRATM(K)=TRATM(K)*DPHINV
     260         END DO
     261         A2=0.0
     262         DO 20 J=I,JN
     263            RR(J)=RM
     264          U(J)=UM
     265          V(J)=VM
     266          DO K=1,NTRA
     267           TRA(J,K)=TRATM(K)
     268          END DO
     269            RDCP=(RD*(1.-RR(J))+RR(J)*RV)/ (CPD*(1.-RR(J))+RR(J)*CPV)
     270            X=(0.001*P(J))**RDCP
     271            T(J)=X
     272            A2=A2+(CPD*(1.-RR(J))+RR(J)*CPV)*X*(PH(J)-PH(J+1))
     273   20    CONTINUE
     274         DO 25 J=I,JN
     275            TH(J)=AHM/A2
     276            T(J)=T(J)*TH(J)
     277   25    CONTINUE
     278   30 CONTINUE
     279
     280      ENDIF ! added sbl
     281C
     282C   ***   RESET INPUT ARRAYS IF ok_adj 0   ***
     283C
     284      IF (ok_adj)THEN
    254285         DO 35 I=1,ND
    255             T1(I)=T(I)
    256             R1(I)=RR(I)
     286
     287           FT2(I)=(T(I)-T1(I))/DELT  ! sbl
     288           FR2(I)=(RR(I)-R1(I))/DELT  ! sbl
     289           FU2(I)=(U(I)-U1(I))/DELT  ! sbl
     290           FV2(I)=(V(I)-V1(I))/DELT  ! sbl
     291
     292c!            T1(I)=T(I)      ! commente sbl
     293c!            R1(I)=RR(I)     ! commente sbl
    257294   35    CONTINUE
    258295      END IF
    … …  
    401438     $        +TV(ICB+1)*(P(ICB)-PBASE   )/(P(ICB)-P(ICB+1))
    402439C
     440c test sb:
     441c@      write(*,*) '++++++++++++++++++++++++++++++++++++++++'
     442c@      write(*,*) 'plcl,dpbas,tvpbas,tvbas,tvp(icb),tvp(icb1)
     443c@     :             ,tv(icb),tv(icb1)'
     444c@      write(*,*) plcl,dpbase,tvpbase,tvbase,tvp(icb)
     445c@     L          ,tvp(icb+1),tv(icb),tv(icb+1)
     446c@      write(*,*) '++++++++++++++++++++++++++++++++++++++++'
     447c fin test sb
    403448      BUOYBASE = TVPBASE - TVBASE
    404449C
    … …  
    418463c sb3d      print *, 'P ',(p(i),i=1,nl)
    419464c sb3d      print *,'ICB ',icb
     465c test sb:
     466c@      write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
     467c@      write(*,*) 'icb,icbs,inb,buoybase:'
     468c@      write(*,*) icb,icb+1,inb,buoybase
     469c@      write(*,*) 'k,tvp,tv,tp,buoy,ep: '
     470c@      do k=1,nl
     471c@      write(*,*) k,tvp(k),tv(k),tp(k),buoy(k),ep(k)
     472c@      enddo
     473c@      write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
     474c fin test sb
     475
     476
    420477C
    421478C   ***   MAKE SURE THAT COLUMN IS DRY ADIABATIC BETWEEN THE SURFACE  ***
    … …  
    552609Cjyg2
    553610c sb3d         print *,'buoy,dlnp,dcape,cape',buoy(i-1),dlnp,dcape,cape
     611c test sb:
     612c@       write(*,*) '############################################'
     613c@         write(*,*) 'cape,rrd,buoy,deltap,p,pbase,ph:'
     614c@     :     ,cape,rd,buoy(i-1),deltap,p(i-1),pbase,ph(i)
     615c@       write(*,*) '############################################'
     616
     617c fin test sb
    554618         CAPE=AMAX1(0.0,CAPE)
    555619C
    … …  
    573637         AMU=0.5*(SIG(I)+SIGOLD)*W
    574638         M(I)=AMU*0.007*P(I)*(PH(I)-PH(I+1))/TV(I)
     639
     640c --------- test sb:
     641c       write(*,*) '############################################'
     642c       write(*,*) 'k,amu,buoy(k-1),deltap,w,beta,fac,cape,w0(k)'
     643c       write(*,*) i,amu,buoy(i-1),deltap
     644c     :           ,w,beta,fac,cape,w0(i)
     645c       write(*,*) '############################################'
     646c ---------
     647
    575648         W0(I)=W
    576649   98 CONTINUE
    … …  
    750823 
    751824***********************************************************
     825c--- test sb:
     826c@       write(*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
     827c@       write(*,*) 'inb,m(inb),ment(inb,inb),sigij(inb,inb):'
     828c@       write(*,*) inb,m(inb),ment(inb,inb),sigij(inb,inb)
     829c@       write(*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
     830c---
     831
    752832 
    753833 
    … …  
    9361016C
    9371017        PRECIP=WT(1)*SIGD*WATER(1)*8640.0
     1018
     1019c sb  ***  Calculate downdraft velocity scale and surface temperature and  ***
     1020c sb  ***                    water vapor fluctuations                      ***
     1021c sb            (inspire de convect 4.3)
     1022
     1023c       BETAD=10.0         
     1024       BETAD=5.0         
     1025       WD=BETAD*ABS(MP(ICB))*0.01*RD*T(ICB)/(SIGD*P(ICB))
     1026
    9381027  405   CONTINUE
    9391028C
    … …  
    10211110         FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I))
    10221111         FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I))
     1112C (saturated updrafts resulting from mixing)      ! cld   
     1113         QCOND(I)=QCOND(I)+(ELIJ(K,I)-AWAT)       ! cld
     1114         NQCOND(I)=NQCOND(I)+1.                   ! cld
    10231115         DO J=1,NTRA
    10241116          FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)-
    … …  
    10451137     1    MP(I)*(TRAP(I,J)-TRAP(I-1,J)))
    10461138        END DO
     1139C (saturated downdrafts resulting from mixing)    ! cld
     1140        DO K=I+1,INB                              ! cld
     1141         QCOND(I)=QCOND(I)+ELIJ(K,I)              ! cld
     1142         NQCOND(I)=NQCOND(I)+1.                   ! cld
     1143        ENDDO                                     ! cld
     1144C (particular case: no detraining level is found) ! cld
     1145        IF (NENT(I).EQ.0) THEN                    ! cld
     1146         QCOND(I)=QCOND(I)+(1-EP(I))*CLW(I)       ! cld
     1147         NQCOND(I)=NQCOND(I)+1.                   ! cld
     1148        ENDIF                                     ! cld
     1149        IF (NQCOND(I).NE.0.) THEN                 ! cld
     1150         QCOND(I)=QCOND(I)/NQCOND(I)              ! cld
     1151        ENDIF                                     ! cld
    10471152  500   CONTINUE
    10481153 
    … …  
    10541159C   ***          INTEGRATED ENTHALPY AND WATER TENDENCIES         ***
    10551160C
     1161c test sb:
     1162c@      write(*,*) '--------------------------------------------'
     1163c@      write(*,*) 'inb,ft,hp,h,t,rr,qent,ment,water,waterp,wt,mp,b'
     1164c@      write(*,*) inb,ft(inb),hp(inb),h(inb)
     1165c@     :   ,t(inb),rr(inb),qent(inb,inb)
     1166c@     :   ,ment(inb,inb),water(inb)
     1167c@     :   ,water(inb+1),wt(inb),mp(inb),b(inb)
     1168c@      write(*,*) '--------------------------------------------'
     1169c fin test sb:
     1170
    10561171      AX=0.1*MENT(INB,INB)*(HP(INB)-H(INB)+T(INB)*
    10571172     1    (CPV-CPD)*(RR(INB)-QENT(INB,INB)))/(CPN(INB)*
    … …  
    11871302       Mke(i)=upwd(i)+dnwd(i)
    11881303       enddo
     1304
     1305C
     1306C   *** Diagnose the in-cloud mixing ratio   ***              ! cld
     1307C   ***           of condensed water         ***              ! cld
     1308C                                                             ! cld
     1309       DO I=1,ND                                              ! cld
     1310        MAA(I)=0.0                                            ! cld
     1311        WA(I)=0.0                                             ! cld
     1312        SIGA(I)=0.0                                           ! cld
     1313       ENDDO                                                  ! cld
     1314       DO I=NK,INB                                            ! cld
     1315       DO K=I+1,INB+1                                         ! cld
     1316        MAA(I)=MAA(I)+M(K)                                    ! cld
     1317       ENDDO                                                  ! cld
     1318       ENDDO                                                  ! cld
     1319       DO I=ICB,INB-1                                         ! cld
     1320        AXC(I)=0.                                             ! cld
     1321        DO J=ICB,I                                            ! cld
     1322         AXC(I)=AXC(I)+RD*(TVP(J)-TV(J))*(PH(J)-PH(J+1))/P(J) ! cld
     1323        ENDDO                                                 ! cld
     1324        IF (AXC(I).GT.0.0) THEN                               ! cld
     1325         WA(I)=SQRT(2.*AXC(I))                                ! cld
     1326        ENDIF                                                 ! cld
     1327       ENDDO                                                  ! cld
     1328       DO I=1,NL                                              ! cld
     1329        IF (WA(I).GT.0.0)                                     ! cld
     1330     1    SIGA(I)=MAA(I)/WA(I)*RD*TVP(I)/P(I)/100./DELTAC     ! cld
     1331        SIGA(I) = MIN(SIGA(I),1.0)                            ! cld
     1332        QCONDC(I)=SIGA(I)*CLW(I)*(1.-EP(I))                   ! cld
     1333     1          + (1.-SIGA(I))*QCOND(I)                       ! cld
     1334       ENDDO                                                  ! cld
     1335
     1336
    11891337c$$$cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    11901338c$$$         call writeg1d(1,klev,ma,'ma  ','ma  ')

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

    r79 r373  
    1       SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,
     1c $Header$
     2c
     3      SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs,
    24     s                   d_t, d_q, d_ql, rneb, radliq, rain, snow,
    3      s                   prfl, psfl)
     5     s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
     6     s                   frac_impa, frac_nucl,
     7     s                   prfl, psfl, rhcl)
     8
    49c
    510      IMPLICIT none
    … …  
    1419#include "dimphy.h"
    1520#include "YOMCST.h"
     21#include "tracstoke.h"
     22#include "fisrtilp.h"
    1623c
    1724c Arguments:
    … …  
    2734      REAL rneb(klon,klev) ! fraction nuageuse
    2835      REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
     36      REAL rhcl(klon,klev) ! humidite relative en ciel clair
    2937      REAL rain(klon) ! pluies (mm/s)
    3038      REAL snow(klon) ! neige (mm/s)
    3139      REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    3240      REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     41cAA
     42c Coeffients de fraction lessivee : pour OFF-LINE
     43c
     44      REAL pfrac_nucl(klon,klev)
     45      REAL pfrac_1nucl(klon,klev)
     46      REAL pfrac_impa(klon,klev)
     47c
     48c Fraction d'aerosols lessivee par impaction et par nucleation
     49c POur ON-LINE
     50c
     51      REAL frac_impa(klon,klev)
     52      REAL frac_nucl(klon,klev)
     53      real zct(klon),zcl(klon)
     54cAA
    3355c
    3456c Options du programme:
    … …  
    3658      REAL seuil_neb ! un nuage existe vraiment au-dela
    3759      PARAMETER (seuil_neb=0.001)
    38       REAL ct ! inverse du temps pour qu'un nuage precipite
    39       PARAMETER (ct=1./1800.)
    40       REAL cl ! seuil de precipitation
    41       PARAMETER (cl=2.6e-4)
    42 ccc      PARAMETER (cl=2.3e-4)
    43 ccc      PARAMETER (cl=2.0e-4)
     60
    4461      INTEGER ninter ! sous-intervals pour la precipitation
    4562      PARAMETER (ninter=5)
    4663      LOGICAL evap_prec ! evaporation de la pluie
    4764      PARAMETER (evap_prec=.TRUE.)
    48       REAL coef_eva
    49       PARAMETER (coef_eva=2.0E-05)
    50       LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
    51       REAL ratqs ! determine la largeur de distribution de vapeur
    52       PARAMETER (calcrat=.TRUE.)
    53       REAL zx_min, rat_max
    54       PARAMETER (zx_min=1.0, rat_max=0.01)
    55       REAL zx_max, rat_min
    56       PARAMETER (zx_max=0.1, rat_min=0.3)
    57       REAL zx
     65      REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
     66      logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
     67
     68      real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
     69      real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
     70      real erf
    5871c
    5972      LOGICAL cpartiel ! condensation partielle
    … …  
    6477c Variables locales:
    6578c
    66       INTEGER i, k, n
     79      INTEGER i, k, n, kk
    6780      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
    6881      REAL zrfl(klon), zrfln(klon), zqev, zqevt
    … …  
    7689      SAVE appel1er
    7790c
     91c---------------------------------------------------------------
     92c
     93cAA Variables traceurs:
     94cAA  Provisoire !!! Parametres alpha du lessivage
     95cAA  A priori on a 4 scavenging # possibles
     96c
     97      REAL a_tr_sca(4)
     98      save a_tr_sca
     99c
     100c Variables intermediaires
     101c
     102      REAL zalpha_tr
     103      REAL zfrac_lessi
     104      REAL zprec_cond(klon)
     105cAA
     106c---------------------------------------------------------------
     107c
    78108c Fonctions en ligne:
    79109c
    80       REAL fallv ! vitesse de chute pour crystaux de glace
     110      REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
    81111      REAL zzz
    82112#include "YOETHF.h"
    83113#include "FCTTRE.h"
    84       fallv (zzz) = 3.29/2.0 * ((zzz)**0.16)
    85 ccc      fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
    86 ccc      fallv (zzz) = 3.29 * ((zzz)**0.16)
     114      fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
     115      fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
    87116c
    88117      DATA appel1er /.TRUE./
    89 c
     118
    90119      IF (appel1er) THEN
    91          PRINT*, 'fisrtilp, calcrat:', calcrat
     120c
    92121         PRINT*, 'fisrtilp, ninter:', ninter
    93122         PRINT*, 'fisrtilp, evap_prec:', evap_prec
    … …  
    96125          PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
    97126          PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
    98           CALL abort
     127c         CALL abort
    99128         ENDIF
    100129         appel1er = .FALSE.
    101       ENDIF
    102 c
     130c
     131cAA initialiation provisoire
     132       a_tr_sca(1) = -0.5
     133       a_tr_sca(2) = -0.5
     134       a_tr_sca(3) = -0.5
     135       a_tr_sca(4) = -0.5
     136c
     137cAA Initialisation a 1 des coefs des fractions lessivees
     138c
     139      DO k = 1, klev
     140       DO i = 1, klon
     141          pfrac_nucl(i,k)=1.
     142          pfrac_1nucl(i,k)=1.
     143          pfrac_impa(i,k)=1.
     144       ENDDO
     145      ENDDO
     146
     147      ENDIF          !  test sur appel1er
     148c
     149cMAf Initialisation a 0 de zoliq
     150       DO i = 1, klon
     151          zoliq(i)=0.
     152       ENDDO
    103153c Determiner les nuages froids par leur temperature
     154c  nexpo regle la raideur de la transition eau liquide / eau glace.
    104155c
    105156      ztglace = RTT - 15.0
    … …  
    115166      ENDDO
    116167      ENDDO
     168
    117169      DO k = 1, klev
    118170      DO i = 1, klon
    … …  
    122174         rneb(i,k) = 0.0
    123175         radliq(i,k) = 0.0
     176         frac_nucl(i,k) = 1.
     177         frac_impa(i,k) = 1.
    124178      ENDDO
    125179      ENDDO
    … …  
    136190      ENDDO
    137191c
     192c
     193cAA Pour plus de securite
     194
     195      zalpha_tr   = 0.
     196      zfrac_lessi = 0.
     197
     198cAA----------------------------------------------------------
     199c
    138200c Boucle verticale (du haut vers le bas)
    139201c
    140202      DO 9999 k = klev, 1, -1
     203c
     204cAA----------------------------------------------------------
    141205c
    142206      DO i = 1, klon
    … …  
    171235         zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
    172236     .                            /RG/dtime
     237
     238c pour la glace, on réévapore toute la précip dans la couche du dessous
     239c la glace venant de la couche du dessus est simplement dans la couche
     240c du dessous.
     241
     242         IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
     243
    173244         zq(i) = zq(i) - (zrfln(i)-zrfl(i))
    174245     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    … …  
    210281c
    211282      IF (cpartiel) THEN
    212          DO i = 1, klon
    213 c
    214             zx = pplay(i,k)/paprs(i,1)
    215             zx = (zx_max-zx)/(zx_max-zx_min)
    216             zx = MIN(MAX(zx,0.0),1.0)
    217             zx = zx * zx * zx
    218             ratqs = zx * (rat_max-rat_min) + rat_min
    219             IF (.NOT.calcrat) ratqs=0.05
    220 ccc            IF (.NOT.calcrat) ratqs=0.2
    221 c
    222             zdelq = ratqs * zq(i)
     283
     284c        print*,'Dans partiel k=',k
     285c
     286c   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
     287c   nuageuse a partir des PDF de Sandrine Bony.
     288c   rneb  : fraction nuageuse
     289c   zqn   : eau totale dans le nuage
     290c   zcond : eau condensee moyenne dans la maille.
     291c           on prend en compte le réchauffement qui diminue la partie condensee
     292c
     293c   Version avec les raqts
     294
     295         if (iflag_pdf.eq.0) then
     296
     297           do i=1,klon
     298            zdelq = min(ratqs(i,k),0.99) * zq(i)
    223299            rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
    224300            zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
     301           enddo
     302
     303         else
     304c
     305c   Version avec les nouvelles PDFs.
     306           do i=1,klon
     307              if(zq(i).lt.1.e-15) then
     308                print*,'ZQ(',i,',',k,')=',zq(i)
     309                zq(i)=1.e-15
     310              endif
     311           enddo
     312           do i=1,klon
     313            zpdf_sig(i)=ratqs(i,k)*zq(i)
     314            zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
     315            zpdf_delta(i)=log(zq(i)/zqs(i))
     316            zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
     317            zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
     318            zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
     319            zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
     320            zpdf_e1(i)=1.-erf(zpdf_e1(i))
     321            zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
     322            zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
     323            zpdf_e2(i)=1.-erf(zpdf_e2(i))
     324            if (zpdf_e1(i).lt.1.e-10) then
     325               rneb(i,k)=0.
     326               zqn(i)=zqs(i)
     327            else
     328               rneb(i,k)=0.5*zpdf_e1(i)
     329               zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
     330            endif
     331           
     332           enddo
     333
     334        endif ! iflag_pdf
     335
     336         do i=1,klon
    225337            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
    226338            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
    227339            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
    228             zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
    229          ENDDO
     340c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
     341c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
     342c  la convection.
     343c  ATTENTION !!! Il va falloir verifier tout ca.
     344            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
     345c           print*,'ZDQS ',zdqs(i)
     346c--Olivier
     347            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
     348            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
     349            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
     350c--fin
     351           ENDDO
    230352      ELSE
    231353         DO i = 1, klon
    … …  
    262384      DO i = 1, klon
    263385      IF (rneb(i,k).GT.0.0) THEN
    264          zchau(i) = ct*dtime/FLOAT(ninter) * zoliq(i)
    265      .          * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) *(1.-zfice(i))
    266386         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
    267          zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
    268      .              *fallv(zrhol(i)) * zfice(i)
     387
     388         if (ptconv(i,k)) then
     389            zcl(i)=cld_lc_con
     390            zct(i)=1./cld_tau_con
     391         else
     392            zcl(i)=cld_lc_lsc
     393            zct(i)=1./cld_tau_lsc
     394         endif
     395c  quantité d'eau à élminier.
     396         zchau(i) = zct(i)*dtime/FLOAT(ninter) * zoliq(i)
     397     .         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl(i))**2)) *(1.-zfice(i))
     398c  meme chose pour la glace.
     399         if (ptconv(i,k)) then
     400            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
     401     .              *fallvc(zrhol(i)) * zfice(i)
     402         else
     403            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
     404     .              *fallvs(zrhol(i)) * zfice(i)
     405         endif
    269406         ztot(i) = zchau(i) + zfroi(i)
    270407         IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
    … …  
    296433      ENDDO
    297434c
     435cAA--------------- Calcul du lessivage stratiforme  -------------
     436
     437      DO i = 1,klon
     438c
     439         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
     440     .                * (paprs(i,k)-paprs(i,k+1))/RG
     441         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
     442cAA lessivage nucleation LMD5 dans la couche elle-meme
     443            if (t(i,k) .GE. ztglace) THEN
     444               zalpha_tr = a_tr_sca(3)
     445            else
     446               zalpha_tr = a_tr_sca(4)
     447            endif
     448            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
     449            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
     450            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
     451c
     452c nucleation avec un facteur -1 au lieu de -0.5
     453            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
     454            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
     455         ENDIF
     456c
     457      ENDDO      ! boucle sur i
     458c
     459cAA Lessivage par impaction dans les couches en-dessous
     460      DO kk = k-1, 1, -1
     461        DO i = 1, klon
     462          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
     463            if (t(i,kk) .GE. ztglace) THEN
     464              zalpha_tr = a_tr_sca(1)
     465            else
     466              zalpha_tr = a_tr_sca(2)
     467            endif
     468            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
     469            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
     470            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
     471          ENDIF
     472        ENDDO
     473      ENDDO
     474c
     475cAA----------------------------------------------------------
     476c                     FIN DE BOUCLE SUR K   
    298477 9999 CONTINUE
     478c
     479cAA-----------------------------------------------------------
    299480c
    300481c Pluie ou neige au sol selon la temperature de la 1ere couche

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

    r353 r373  
    77     .           fder,radsol,frugs,agesno,clesphy0,
    88     .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,tabcntr0,
    9      .           t_ancien,q_ancien,ancien_ok)
     9     .           t_ancien,q_ancien,ancien_ok, rnebcon, ratqs,clwcon)
    1010      IMPLICIT none
    1111c======================================================================
    … …  
    5454
    5555      REAL t_ancien(klon,klev), q_ancien(klon,klev)
     56      real rnebcon(klon,klev),clwcon(klon,klev),ratqs(klon,klev)
    5657      LOGICAL ancien_ok
    5758
    … …  
    11701171      ENDIF
    11711172c
     1173      ierr = NF_INQ_VARID (nid, "CLWCON", nvarid)
     1174      IF (ierr.NE.NF_NOERR) THEN
     1175         PRINT*, "phyetat0: Le champ CLWCON est absent"
     1176         PRINT*, "Depart legerement fausse. Mais je continue"
     1177         clwcon = 0.
     1178      ELSE
     1179#ifdef NC_DOUBLE
     1180         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, clwcon)
     1181#else
     1182         ierr = NF_GET_VAR_REAL(nid, nvarid, clwcon)
     1183#endif
     1184         IF (ierr.NE.NF_NOERR) THEN
     1185            PRINT*, "phyetat0: Lecture echouee pour <CLWCON>"
     1186            CALL abort
     1187         ENDIF
     1188      ENDIF
     1189      xmin = 1.0E+20
     1190      xmax = -1.0E+20
     1191      xmin = MINval(clwcon)
     1192      xmax = MAXval(clwcon)
     1193      PRINT*,'Eau liquide convective (ecart-type) clwcon:', xmin, xmax
     1194c
     1195      ierr = NF_INQ_VARID (nid, "RNEBCON", nvarid)
     1196      IF (ierr.NE.NF_NOERR) THEN
     1197         PRINT*, "phyetat0: Le champ RNEBCON est absent"
     1198         PRINT*, "Depart legerement fausse. Mais je continue"
     1199         rnebcon = 0.
     1200      ELSE
     1201#ifdef NC_DOUBLE
     1202         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rnebcon)
     1203#else
     1204         ierr = NF_GET_VAR_REAL(nid, nvarid, rnebcon)
     1205#endif
     1206         IF (ierr.NE.NF_NOERR) THEN
     1207            PRINT*, "phyetat0: Lecture echouee pour <RNEBCON>"
     1208            CALL abort
     1209         ENDIF
     1210      ENDIF
     1211      xmin = 1.0E+20
     1212      xmax = -1.0E+20
     1213      xmin = MINval(rnebcon)
     1214      xmax = MAXval(rnebcon)
     1215      PRINT*,'Nebulosite convective (ecart-type) rnebcon:', xmin, xmax
     1216
     1217c
     1218      ierr = NF_INQ_VARID (nid, "QANCIEN", nvarid)
     1219      IF (ierr.NE.NF_NOERR) THEN
     1220         PRINT*, "phyetat0: Le champ <QANCIEN> est absent"
     1221         PRINT*, "Depart legerement fausse. Mais je continue"
     1222         ancien_ok = .FALSE.
     1223      ELSE
     1224#ifdef NC_DOUBLE
     1225         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, q_ancien)
     1226#else
     1227         ierr = NF_GET_VAR_REAL(nid, nvarid, q_ancien)
     1228#endif
     1229         IF (ierr.NE.NF_NOERR) THEN
     1230            PRINT*, "phyetat0: Lecture echouee pour <QANCIEN>"
     1231            CALL abort
     1232         ENDIF
     1233      ENDIF
     1234c
     1235      ierr = NF_INQ_VARID (nid, "RATQS", nvarid)
     1236      IF (ierr.NE.NF_NOERR) THEN
     1237         PRINT*, "phyetat0: Le champ <RATQS> est absent"
     1238         PRINT*, "Depart legerement fausse. Mais je continue"
     1239         ratqs = 0.
     1240      ELSE
     1241#ifdef NC_DOUBLE
     1242         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ratqs)
     1243#else
     1244         ierr = NF_GET_VAR_REAL(nid, nvarid, ratqs)
     1245#endif
     1246         IF (ierr.NE.NF_NOERR) THEN
     1247            PRINT*, "phyetat0: Lecture echouee pour <RATQS>"
     1248            CALL abort
     1249         ENDIF
     1250      ENDIF
     1251      xmin = 1.0E+20
     1252      xmax = -1.0E+20
     1253      xmin = MINval(ratqs)
     1254      xmax = MAXval(ratqs)
     1255      PRINT*,'(ecart-type) ratqs:', xmin, xmax
     1256c
    11721257c Fermer le fichier:
    11731258c

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

    r352 r373  
    77     .           radsol,frugs,agesno,
    88     .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,
    9      .           t_ancien, q_ancien)
     9     .           t_ancien, q_ancien, rnebcon, ratqs, clwcon)
    1010      IMPLICIT none
    1111c======================================================================
    … …  
    5353      REAL pctsrf(klon, nbsrf)
    5454      REAL t_ancien(klon,klev), q_ancien(klon,klev)
     55      real clwcon(klon,klev),rnebcon(klon,klev),ratqs(klon,klev)
    5556c
    5657      INTEGER nid, nvarid, idim1, idim2, idim3
    … …  
    9697      IF(     ok_orolf ) tab_cntrl(11 ) = 1.
    9798
    98       tab_cntrl(13) = dayref
    99       tab_cntrl(14) = anneeref
    10099      tab_cntrl(13) = day_end
    101100      tab_cntrl(14) = annee_ref
    … …  
    670669#endif
    671670c
     671      ierr = NF_REDEF (nid)
     672#ifdef NC_DOUBLE
     673      ierr = NF_DEF_VAR (nid, "CLWCON", NF_DOUBLE, 1, idim2,nvarid)
     674#else
     675      ierr = NF_DEF_VAR (nid, "CLWCON", NF_FLOAT, 1, idim2,nvarid)
     676#endif
     677      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
     678     .                        "Eau liquide convective")
     679      ierr = NF_ENDDEF(nid)
     680#ifdef NC_DOUBLE
     681      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,clwcon)
     682#else
     683      ierr = NF_PUT_VAR_REAL (nid,nvarid,clwcon)
     684#endif
     685c
     686      ierr = NF_REDEF (nid)
     687#ifdef NC_DOUBLE
     688      ierr = NF_DEF_VAR (nid, "RNEBCON", NF_DOUBLE, 1, idim2,nvarid)
     689#else
     690      ierr = NF_DEF_VAR (nid, "RNEBCON", NF_FLOAT, 1, idim2,nvarid)
     691#endif
     692      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
     693     .                        "Nebulosite convective")
     694      ierr = NF_ENDDEF(nid)
     695#ifdef NC_DOUBLE
     696      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rnebcon)
     697#else
     698      ierr = NF_PUT_VAR_REAL (nid,nvarid,rnebcon)
     699#endif
     700c
     701      ierr = NF_REDEF (nid)
     702#ifdef NC_DOUBLE
     703      ierr = NF_DEF_VAR (nid, "RATQS", NF_DOUBLE, 1, idim2,nvarid)
     704#else
     705      ierr = NF_DEF_VAR (nid, "RATQS", NF_FLOAT, 1, idim2,nvarid)
     706#endif
     707      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
     708     .                        "Ratqs")
     709      ierr = NF_ENDDEF(nid)
     710#ifdef NC_DOUBLE
     711      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ratqs)
     712#else
     713      ierr = NF_PUT_VAR_REAL (nid,nvarid,ratqs)
     714#endif
     715c
    672716c
    673717      ierr = NF_CLOSE(nid)

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

    r364 r373  
    388388      EXTERNAL conlmd    ! convection (schema LMD)
    389389cKE43
    390       EXTERNAL conema  ! convect4.3
     390      EXTERNAL conema3  ! convect4.3
    391391      EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
    392392cAA
    … …  
    412412c Variables locales
    413413c
     414      real clwcon(klon,klev),rnebcon(klon,klev)
     415      real clwcon0(klon,klev),rnebcon0(klon,klev)
     416      save rnebcon, clwcon
     417
     418      REAL rhcl(klon,klev)    ! humiditi relative ciel clair
    414419      REAL dialiq(klon,klev)  ! eau liquide nuageuse
    415420      REAL diafra(klon,klev)  ! fraction nuageuse
    … …  
    464469c
    465470      REAL za, zb
    466       REAL zx_t, zx_qs, zdelta, zcor, zlvdcp, zlsdcp
    467       INTEGER i, k, iq, nsrf, ll
     471      REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
     472      real zqsat(klon,klev)
     473      INTEGER i, k, iq, ig, j, nsrf, ll
    468474      REAL t_coup
    469475      PARAMETER (t_coup=234.0)
    … …  
    534540      REAL d_t_lif(klon,klev)
    535541
    536       REAL ratqs(klon,klev)
    537       integer flag_ratqs
     542      REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
     543      real ratqsbas,ratqshaut
     544      save ratqsbas,ratqshaut, ratqs
    538545      real zpt_conv(klon,klev)
     546
     547c Parametres lies au nouveau schema de nuages (SB, PDF)
     548      real fact_cldcon
     549      real facttemps
     550      logical ok_newmicro
     551      save ok_newmicro
     552      save fact_cldcon,facttemps
     553
     554      integer iflag_cldcon
     555      save iflag_cldcon
     556
     557      logical ptconv(klon,klev)
    539558
    540559c
    … …  
    630649c
    631650         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
    632      .                  ok_instan)
     651     .                  ok_instan, fact_cldcon, facttemps,ok_newmicro,
     652     .                  iflag_cldcon,ratqsbas,ratqshaut)
    633653
    634654         DO k = 2, nvm          ! pas de vegetation
    … …  
    655675     .       dlw,radsol,frugs,agesno,clesphy0,
    656676     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
    657      .       t_ancien, q_ancien, ancien_ok )
     677     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon )
    658678
    659679c
    … …  
    13071327     .                "ave(X)", zsto,zout)
    13081328c
     1329         CALL histdef(nid_mth, "ducon", "Convection du", "m/s2",
     1330     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     1331     .                "ave(X)", zsto,zout)
     1332c
    13091333         CALL histdef(nid_mth, "dqcon", "Convection dQ", "Kg/Kg/s",
    13101334     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
    … …  
    17741798      DO i = 1, klon
    17751799         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
    1776          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
     1800c        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
     1801         zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
    17771802         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
    17781803         zb = MAX(0.0,ql_seri(i,k))
    … …  
    19812006          else
    19822007
    1983           CALL conema (dtime,paprs,pplay,t_seri,q_seri,
     2008c          print*,'Avant conema OUI'
     2009          CALL conema3 (dtime,
     2010     .        paprs,pplay,t_seri,q_seri,
    19842011     .        u_seri,v_seri,tr_seri,nbtr,
    19852012     .        ema_work1,ema_work2,
    … …  
    19882015     .        upwd,dnwd,dnwd0,bas,top,
    19892016     .        Ma,cape,tvp,rflag,
    1990      .       pbase
    1991      .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
     2017     .        pbase
     2018     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
     2019     .        ,clwcon0)
     2020c          print*,'Apres conema3 '
     2021
     2022c Calculer l'humidite relative pour diagnostique
     2023c
     2024      DO k = 1, klev
     2025      DO i = 1, klon
     2026         zx_t = t_seri(i,k)
     2027         IF (thermcep) THEN
     2028            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
     2029            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
     2030            zx_qs  = MIN(0.5,zx_qs)
     2031            zcor   = 1./(1.-retv*zx_qs)
     2032            zx_qs  = zx_qs*zcor
     2033         ELSE
     2034           IF (zx_t.LT.t_coup) THEN
     2035              zx_qs = qsats(zx_t)/pplay(i,k)
     2036           ELSE
     2037              zx_qs = qsatl(zx_t)/pplay(i,k)
     2038           ENDIF
     2039         ENDIF
     2040         zqsat(i,k)=zx_qs
     2041      ENDDO
     2042      ENDDO
     2043
     2044c   calcul des propriétés des nuages convectifs
     2045             clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
     2046             call clouds_gno
     2047     s       (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
     2048
    19922049          endif
    19932050          DO i = 1, klon
    … …  
    20542111      IF (nqmax.GT.2) THEN !--melange convectif de traceurs
    20552112c
    2056           IF (iflag_con .LT. 2 .AND.  iflag_con .GT. 4 ) THEN
     2113          IF (iflag_con .NE. 2 .AND. debut) THEN
    20572114              PRINT*, 'Pour l instant, seul conflx fonctionne ',
    20582115     $            'avec traceurs', iflag_con
    20592116              PRINT*,' Mettre iflag_con',
    2060      $            ' = 2, 3 ou 4 dans run.def et repasser'
    2061               CALL abort
     2117     $            ' = 2 dans run.def et repasser'
     2118c              CALL abort
    20622119              ENDIF
    20632120c
    … …  
    20742131      ENDDO
    20752132
    2076 c   RATQS
    2077       if (iflag_con.eq.2) then
    2078           flag_ratqs=0
     2133
     2134c-------------------------------------------------------------------------
     2135c  Caclul des ratqs
     2136c-------------------------------------------------------------------------
     2137
     2138c      print*,'calcul des ratqs'
     2139c   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
     2140c   ----------------
     2141c   on ecrase le tableau ratqsc calcule par clouds_gno
     2142      if (iflag_cldcon.eq.1) then
     2143         do k=1,klev
     2144         do i=1,klon
     2145            if(ptconv(i,k)) then
     2146              ratqsc(i,k)=ratqsbas
     2147     s        +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
     2148            else
     2149               ratqsc(i,k)=0.
     2150            endif
     2151         enddo
     2152         enddo
     2153      endif
     2154
     2155c   ratqs stables
     2156c   -------------
     2157      do k=1,klev
     2158         ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
     2159     s   min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)       
     2160      enddo
     2161
     2162
     2163c  ratqs final
     2164c  -----------
     2165      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
     2166c   les ratqs sont une conbinaison de ratqss et ratqsc
     2167c   ratqs final
     2168c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
     2169c   relaxation des ratqs
     2170         facttemps=exp(-pdtphys/1.e4)
     2171         ratqs(:,:)=max(ratqs(:,:)*facttemps,ratqss(:,:))
     2172         ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
     2173c         print*,'calcul des ratqs fini'
    20792174      else
    2080           flag_ratqs=1
     2175c   on ne prend que le ratqs stable pour fisrtilp
     2176         ratqs(:,:)=ratqss(:,:)
    20812177      endif
    2082       call calcratqs (flag_ratqs,
    2083      I            paprs,pplay,q_seri,d_t_con,d_t_ajs
    2084      O           ,ratqs,zpt_conv)
     2178
     2179
    20852180c
    20862181c Appeler le processus de condensation a grande echelle
    20872182c et le processus de precipitation
    2088 c
    2089       CALL fisrtilp_tr(dtime,paprs,pplay,
    2090      .           t_seri, q_seri,ratqs,
     2183c-------------------------------------------------------------------------
     2184      CALL fisrtilp(dtime,paprs,pplay,
     2185     .           t_seri, q_seri,ptconv,ratqs,
    20912186     .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
    20922187     .           rain_lsc, snow_lsc,
    20932188     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
    20942189     .           frac_impa, frac_nucl,
    2095      .           prfl, psfl)
     2190     .           prfl, psfl, rhcl)
     2191
    20962192      WHERE (rain_lsc < 0) rain_lsc = 0.
    20972193      WHERE (snow_lsc < 0) snow_lsc = 0.
    … …  
    21182214      ENDIF
    21192215c
    2120 c Nuages diagnostiques:
    2121 c
    2122       IF (iflag_con.EQ.2) THEN ! seulement pour Tiedtke
     2216c-------------------------------------------------------------------
     2217c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
     2218c-------------------------------------------------------------------
     2219
     2220c 1. NUAGES CONVECTIFS
     2221c
     2222      IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
     2223
     2224c Nuages diagnostiques pour Tiedtke
    21232225      CALL diagcld1(paprs,pplay,
    21242226     .             rain_con,snow_con,ibas_con,itop_con,
    … …  
    21322234      ENDDO
    21332235      ENDDO
     2236
     2237      ELSE IF (iflag_cldcon.eq.3) THEN
     2238c  On prend pour les nuages convectifs le max du calcul de la
     2239c  convection et du calcul du pas de temps précédent diminué d'un facteur
     2240c  facttemps
     2241      facttemps=pdtphys/1.e4
     2242      do k=1,klev
     2243         do i=1,klon
     2244            rnebcon(i,k)=rnebcon(i,k)*facttemps
     2245            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
     2246     s      then
     2247                rnebcon(i,k)=rnebcon0(i,k)
     2248                clwcon(i,k)=clwcon0(i,k)
     2249            endif
     2250         enddo
     2251      enddo
     2252
     2253c   On prend la somme des fractions nuageuses et des contenus en eau
     2254      cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
     2255      cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
     2256
     2257
    21342258      ENDIF
    21352259c
    2136 c Nuages stratus artificiels:
     2260c 2. NUAGES STARTIFORMES
    21372261c
    21382262      IF (ok_stratus) THEN
    … …  
    21742298         ENDIF
    21752299         zx_rh(i,k) = q_seri(i,k)/zx_qs
     2300         zqsat(i,k)=zx_qs
    21762301      ENDDO
    21772302      ENDDO
    … …  
    21802305c parametres pour diagnostiques:
    21812306c
     2307      if (ok_newmicro) then
     2308      CALL newmicro (paprs, pplay,ok_newmicro,
     2309     .            t_seri, cldliq, cldfra, cldtau, cldemi,
     2310     .            cldh, cldl, cldm, cldt, cldq)
     2311      else
    21822312      CALL nuage (paprs, pplay,
    21832313     .            t_seri, cldliq, cldfra, cldtau, cldemi,
    21842314     .            cldh, cldl, cldm, cldt, cldq)
     2315      endif
    21852316c
    21862317c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
    … …  
    33613492     .      radsol,frugs,agesno,
    33623493     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
    3363      .      t_ancien, q_ancien)
     3494     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon)
    33643495      ENDIF
    33653496
Note: See TracChangeset for help on using the changeset viewer.