Ignore:
Timestamp:
Jul 9, 2014, 11:19:07 PM (10 years ago)
Author:
fhourdin
Message:

nclusion de la thermodynamique de la glace
Ice thermodynamics
(Catherine Rio)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/trunk/libf/phylmd/physiq.F90

    r2076 r2086  
    201201  INTEGER iliq          ! indice de traceurs pour eau liquide
    202202  PARAMETER (iliq=2)
    203 
     203!CR: on ajoute la phase glace
     204  INTEGER isol          ! indice de traceurs pour eau glace
     205  PARAMETER (isol=3)
    204206  !
    205207  !
     
    596598
    597599  ! tendance nulles
    598   REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0
     600  REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
    599601
    600602  !
     
    13771379  dq0(:,:)=0.
    13781380  dql0(:,:)=0.
     1381  dqi0(:,:)=0.
    13791382  !
    13801383  ! Mettre a zero des variables de sortie (pour securite)
     
    14221425        q_seri(i,k)  = qx(i,k,ivap)
    14231426        ql_seri(i,k) = qx(i,k,iliq)
    1424         qs_seri(i,k) = 0.
     1427!CR: ATTENTION, on rajoute la variable glace
     1428        if (nqo.eq.2) then
     1429           qs_seri(i,k) = 0.
     1430        else if (nqo.eq.3) then
     1431           qs_seri(i,k) = qx(i,k,isol)
     1432        endif
    14251433     ENDDO
    14261434  ENDDO
    14271435  tke0(:,:)=pbl_tke(:,:,is_ave)
    1428   IF (nqtot.GE.3) THEN
    1429      DO iq = 3, nqtot
     1436!CR:Nombre de traceurs de l'eau: nqo
     1437!  IF (nqtot.GE.3) THEN
     1438   IF (nqtot.GE.(nqo+1)) THEN
     1439!     DO iq = 3, nqtot       
     1440     DO iq = nqo+1, nqtot 
    14301441        DO  k = 1, klev
    14311442           DO  i = 1, klon
    1432               tr_seri(i,k,iq-2) = qx(i,k,iq)
     1443!              tr_seri(i,k,iq-2) = qx(i,k,iq)
     1444              tr_seri(i,k,iq-nqo) = qx(i,k,iq)
    14331445           ENDDO
    14341446        ENDDO
     
    16171629        ENDIF
    16181630        !>jyg
     1631     
     1632        if (iflag_ice_thermo.eq.0) then   
     1633!pas necessaire a priori
    16191634
    16201635        zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
     
    16271642        d_t_eva(i,k) = za
    16281643        d_q_eva(i,k) = zb
     1644
     1645        else
     1646
     1647!CR: on ré-évapore eau liquide et glace
     1648
     1649!        zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
     1650!        zb = MAX(0.0,ql_seri(i,k))
     1651!        za = - MAX(0.0,ql_seri(i,k)) &
     1652!             * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
     1653        zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
     1654        za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
     1655             - MAX(0.0,qs_seri(i,k))*zlsdcp
     1656        t_seri(i,k) = t_seri(i,k) + za
     1657        q_seri(i,k) = q_seri(i,k) + zb
     1658        ql_seri(i,k) = 0.0
     1659!on évapore la glace
     1660        qs_seri(i,k) = 0.0
     1661        d_t_eva(i,k) = za
     1662        d_q_eva(i,k) = zb
     1663        endif
     1664
    16291665     ENDDO
    16301666  ENDDO
     
    17641800     IF (klon_glo==1) THEN
    17651801        CALL add_pbl_tend &
    1766              (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
     1802             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
    17671803     ELSE
    17681804        CALL add_phys_tend &
    1769              (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
     1805             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
    17701806     ENDIF
    17711807     !--------------------------------------------------------------------
     
    21352171  !-----------------------------------------------------------------------------------------
    21362172  ! ajout des tendances de la diffusion turbulente
    2137   CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,paprs,'con')
     2173  CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,dqi0,paprs,'con')
    21382174  !-----------------------------------------------------------------------------------------
    21392175
     
    22522288     d_t_wake(:,:)=dt_wake(:,:)*dtime
    22532289     d_q_wake(:,:)=dq_wake(:,:)*dtime
    2254      CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,paprs,'wake')
     2290     CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake')
    22552291     !-----------------------------------------------------------------------------------------
    22562292
     
    25122548        !-----------------------------------------------------------------------------------------
    25132549        ! ajout des tendances de l'ajustement sec ou des thermiques
    2514         CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,paprs,'ajsb')
     2550        CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs,'ajsb')
    25152551        d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
    25162552        d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
     
    25582594  CALL fisrtilp(dtime,paprs,pplay, &
    25592595       t_seri, q_seri,ptconv,ratqs, &
    2560        d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
     2596       d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
    25612597       rain_lsc, snow_lsc, &
    25622598       pfrac_impa, pfrac_nucl, pfrac_1nucl, &
     
    25702606  !-----------------------------------------------------------------------------------------
    25712607  ! ajout des tendances de la diffusion turbulente
    2572   CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,paprs,'lsc')
     2608  CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs,'lsc')
    25732609  !-----------------------------------------------------------------------------------------
    25742610  DO k = 1, klev
    25752611     DO i = 1, klon
    25762612        cldfra(i,k) = rneb(i,k)
     2613!CR: a quoi ca sert? Faut-il ajouter qs_seri?
    25772614        IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
    25782615     ENDDO
     
    28452882        zx_t = t_seri(i,k)
    28462883        IF (thermcep) THEN
     2884           if (iflag_ice_thermo.eq.0) then
    28472885           zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
     2886           else
     2887           zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))
     2888           endif
    28482889           zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
    28492890           zx_qs  = MIN(0.5,zx_qs)
     
    32883329     !-----------------------------------------------------------------------------------------
    32893330     ! ajout des tendances de la trainee de l'orographie
    3290      CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,paprs,'oro')
     3331     CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro')
    32913332     !-----------------------------------------------------------------------------------------
    32923333     !
     
    33343375     !-----------------------------------------------------------------------------------------
    33353376     ! ajout des tendances de la portance de l'orographie
    3336      CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,paprs,'lif')
     3377     CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,dqi0,paprs,'lif')
    33373378     !-----------------------------------------------------------------------------------------
    33383379     !
     
    33483389     !
    33493390     !  ajout des tendances
    3350      CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,paprs,'hin')
     3391     CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,dqi0,paprs,'hin')
    33513392
    33523393  ENDIF
     
    33563397          rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
    33573398          du_gwd_rando, dv_gwd_rando)
    3358      CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,paprs, &
     3399     CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,dqi0,paprs, &
    33593400          'flott_gwd_rando')
    33603401  end if
     
    36233664        d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
    36243665        d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
     3666!CR: on ajoute le contenu en glace
     3667        if (nqo.eq.3) then
     3668        d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
     3669        endif
    36253670     ENDDO
    36263671  ENDDO
    36273672  !
    3628   IF (nqtot.GE.3) THEN
    3629      DO iq = 3, nqtot
     3673!CR: nb de traceurs eau: nqo
     3674!  IF (nqtot.GE.3) THEN
     3675   IF (nqtot.GE.(nqo+1)) THEN
     3676!     DO iq = 3, nqtot
     3677     DO iq = nqo+1, nqtot
    36303678        DO  k = 1, klev
    36313679           DO  i = 1, klon
    3632               d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
     3680!              d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
     3681               d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
    36333682           ENDDO
    36343683        ENDDO
     
    36533702
    36543703!!! RomP >>>
    3655   IF (nqtot.GE.3) THEN
    3656      DO iq = 3, nqtot
     3704!CR: nb de traceurs eau: nqo
     3705!  IF (nqtot.GE.3) THEN
     3706   IF (nqtot.GE.(nqo+1)) THEN
     3707!     DO iq = 3, nqtot
     3708     DO iq = nqo+1, nqtot
    36573709        DO k = 1, klev
    36583710           DO i = 1, klon
    3659               tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
     3711!              tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
     3712              tr_ancien(i,k,iq-nqo) = tr_seri(i,k,iq-nqo)
    36603713           ENDDO
    36613714        ENDDO
Note: See TracChangeset for help on using the changeset viewer.