Ignore:
Timestamp:
Aug 3, 2024, 2:56:58 PM (4 months ago)
Author:
abarral
Message:

Put .h into modules

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/lsc_scav_orig.F90

    r5159 r5160  
    11!$Id $
    22
    3 SUBROUTINE lsc_scav_orig(pdtime,it,iflag_lscav,oliq,flxr,flxs,rneb,beta_fisrt, &
    4                     beta_v1,pplay,paprs,t,tr_seri,d_tr_insc,          &
    5                     d_tr_bcscav,d_tr_evap,qPrls)
     3SUBROUTINE lsc_scav_orig(pdtime, it, iflag_lscav, oliq, flxr, flxs, rneb, beta_fisrt, &
     4        beta_v1, pplay, paprs, t, tr_seri, d_tr_insc, &
     5        d_tr_bcscav, d_tr_evap, qPrls)
    66  USE ioipsl
    77  USE dimphy
    … …  
    99  USE lmdz_phys_para
    1010  USE traclmdz_mod
    11   USE infotrac,ONLY: nbtr
     11  USE infotrac, ONLY: nbtr
    1212  USE iophy
    1313  USE lmdz_yomcst
    1414  USE lmdz_YOECUMF
    15 
    1615  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
     16  USE lmdz_chem, ONLY: idms, iso2, iso4, ih2s, idmso, imsa, ih2o2, &
     17          n_avogadro, masse_s, masse_so4, rho_water, rho_ice
     18
    1719  IMPLICIT NONE
    18 !=====================================================================
    19 ! Objet : depot humide (lessivage et evaporation) de traceurs
    20 ! Inspired by routines of Olivier Boucher (mars 1998)
    21 ! author R. Pilon 10 octobre 2012
    22 ! last modification 16/01/2013 (reformulation partie evaporation)
    23 !=====================================================================
    24 
    25 
    26   INCLUDE "chem.h"
    27 
    28   REAL,INTENT(IN)                        :: pdtime ! time step (s)
    29   INTEGER,INTENT(IN)                     :: it     ! tracer number
    30   INTEGER,INTENT(IN)                     :: iflag_lscav ! LS scavenging param:
    31 !                                             3=Reddy_Boucher2004, 4=3+RPilon.
    32   REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxr     ! flux precipitant de pluie
    33   REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxs     ! flux precipitant de neige
    34   REAL,INTENT(IN)                        :: oliq ! contenu en eau liquide dans le nuage (kg/kg)
    35   REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb
    36   REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay    ! pression
    37   REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs    ! pression
    38   REAL,DIMENSION(klon,klev),INTENT(IN)   :: t        ! temperature
    39 ! tracers
    40   REAL,DIMENSION(klon,klev,nbtr),INTENT(IN)   :: tr_seri        ! q de traceur 
    41   REAL,DIMENSION(klon,klev),INTENT(IN)        :: beta_fisrt     ! taux de conversion de l'eau cond
    42   REAL,DIMENSION(klon,klev),INTENT(OUT)       :: beta_v1        ! -- (originale version)
    43   REAL,DIMENSION(klon)                        :: his_dh         ! tendance de traceur integre verticalement
    44   REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)  :: d_tr_insc      ! tendance du traceur
    45   REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)  :: d_tr_bcscav  ! tendance de traceur
    46   REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)  :: d_tr_evap
    47   REAL,DIMENSION(klon,nbtr),INTENT(OUT)       :: qPrls      !jyg: concentration tra dans pluie LS a la surf.
    48   REAL :: dxin,dxev                              ! tendance temporaire de traceur incloud
    49   REAL,DIMENSION(klon,klev) :: dxbc       ! tendance temporaire de traceur bc
    50 
    51 
    52 !  variables locales     
    53  LOGICAL,SAVE :: debut=.TRUE.
    54 !$OMP THREADPRIVATE(debut)
    55 
    56   REAL,PARAMETER :: henry=1.4  ! constante de Henry en mol/l/atm ~1.4 for gases
    57   REAL           :: henry_t    !  constante de Henry a T t  (mol/l/atm)
    58   REAL,PARAMETER :: kk=2900.   ! coefficient de dependence en T (K)
     20  !=====================================================================
     21  ! Objet : depot humide (lessivage et evaporation) de traceurs
     22  ! Inspired by routines of Olivier Boucher (mars 1998)
     23  ! author R. Pilon 10 octobre 2012
     24  ! last modification 16/01/2013 (reformulation partie evaporation)
     25  !=====================================================================
     26
     27  REAL, INTENT(IN) :: pdtime ! time step (s)
     28  INTEGER, INTENT(IN) :: it     ! tracer number
     29  INTEGER, INTENT(IN) :: iflag_lscav ! LS scavenging param:
     30  !                                             3=Reddy_Boucher2004, 4=3+RPilon.
     31  REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: flxr     ! flux precipitant de pluie
     32  REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: flxs     ! flux precipitant de neige
     33  REAL, INTENT(IN) :: oliq ! contenu en eau liquide dans le nuage (kg/kg)
     34  REAL, DIMENSION(klon, klev), INTENT(IN) :: rneb
     35  REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay    ! pression
     36  REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs    ! pression
     37  REAL, DIMENSION(klon, klev), INTENT(IN) :: t        ! temperature
     38  ! tracers
     39  REAL, DIMENSION(klon, klev, nbtr), INTENT(IN) :: tr_seri        ! q de traceur
     40  REAL, DIMENSION(klon, klev), INTENT(IN) :: beta_fisrt     ! taux de conversion de l'eau cond
     41  REAL, DIMENSION(klon, klev), INTENT(OUT) :: beta_v1        ! -- (originale version)
     42  REAL, DIMENSION(klon) :: his_dh         ! tendance de traceur integre verticalement
     43  REAL, DIMENSION(klon, klev, nbtr), INTENT(OUT) :: d_tr_insc      ! tendance du traceur
     44  REAL, DIMENSION(klon, klev, nbtr), INTENT(OUT) :: d_tr_bcscav  ! tendance de traceur
     45  REAL, DIMENSION(klon, klev, nbtr), INTENT(OUT) :: d_tr_evap
     46  REAL, DIMENSION(klon, nbtr), INTENT(OUT) :: qPrls      !jyg: concentration tra dans pluie LS a la surf.
     47  REAL :: dxin, dxev                              ! tendance temporaire de traceur incloud
     48  REAL, DIMENSION(klon, klev) :: dxbc       ! tendance temporaire de traceur bc
     49
     50
     51  !  variables locales
     52  LOGICAL, SAVE :: debut = .TRUE.
     53  !$OMP THREADPRIVATE(debut)
     54
     55  REAL, PARAMETER :: henry = 1.4  ! constante de Henry en mol/l/atm ~1.4 for gases
     56  REAL :: henry_t    !  constante de Henry a T t  (mol/l/atm)
     57  REAL, PARAMETER :: kk = 2900.   ! coefficient de dependence en T (K)
    5958  REAL :: f_a     !  rapport de la phase aqueuse a la phase gazeuse
    6059  REAL :: beta    !  taux de conversion de l'eau en pluie
    6160
    6261  INTEGER :: i, k
    63   REAL,DIMENSION(klon,klev)    :: scav  !  water liquid content / fraction aqueuse du constituant
    64   REAL,DIMENSION(klon,klev)    :: zrho
    65   REAL,DIMENSION(klon,klev)    :: zdz
    66   REAL,DIMENSION(klon,klev)    :: zmass ! layer mass
    67 
    68   REAL           :: frac_ev       ! cste pour la reevaporation : dropplet shrinking
    69 !  frac_ev = frac_gas ou frac_aer
    70   REAL,PARAMETER :: frac_gas=1.0  ! cste pour la reevaporation pour les gaz
    71   REAL           :: frac_aer      ! cste pour la reevaporation pour les particules
    72   REAL,DIMENSION(klon,klev) :: deltaP     ! P(i+1)-P(i)
    73   REAL,DIMENSION(klon,klev) :: beta_ev    !  dP/P(i+1)
    74 
    75 !  101.325  m3/l x Pa/atm
    76 !  R        Pa.m3/mol/K
    77 !   cste de dissolution pour le depot humide
    78   REAL,SAVE :: frac_fine_scav
    79   REAL,SAVE :: frac_coar_scav
    80 !$OMP THREADPRIVATE(frac_fine_scav, frac_coar_scav)
    81 
    82 ! below-cloud scav variables
    83 ! aerosol : alpha_r=0.001, gas 0.001  (Pruppacher & Klett 1967)
    84   REAL,SAVE :: alpha_r  !  coefficient d'impaction pour la pluie
    85   REAL,SAVE :: alpha_s  !  coefficient d'impaction pour la neige 
    86   REAL,SAVE :: R_r      !  mean raindrop radius (m)
    87   REAL,SAVE :: R_s      !  mean snow crystal radius (m)
    88 !$OMP THREADPRIVATE(alpha_r, alpha_s, R_r, R_s)
    89   REAL           :: pr, ps, ice, water
     62  REAL, DIMENSION(klon, klev) :: scav  !  water liquid content / fraction aqueuse du constituant
     63  REAL, DIMENSION(klon, klev) :: zrho
     64  REAL, DIMENSION(klon, klev) :: zdz
     65  REAL, DIMENSION(klon, klev) :: zmass ! layer mass
     66
     67  REAL :: frac_ev       ! cste pour la reevaporation : dropplet shrinking
     68  !  frac_ev = frac_gas ou frac_aer
     69  REAL, PARAMETER :: frac_gas = 1.0  ! cste pour la reevaporation pour les gaz
     70  REAL :: frac_aer      ! cste pour la reevaporation pour les particules
     71  REAL, DIMENSION(klon, klev) :: deltaP     ! P(i+1)-P(i)
     72  REAL, DIMENSION(klon, klev) :: beta_ev    !  dP/P(i+1)
     73
     74  !  101.325  m3/l x Pa/atm
     75  !  R        Pa.m3/mol/K
     76  !   cste de dissolution pour le depot humide
     77  REAL, SAVE :: frac_fine_scav
     78  REAL, SAVE :: frac_coar_scav
     79  !$OMP THREADPRIVATE(frac_fine_scav, frac_coar_scav)
     80
     81  ! below-cloud scav variables
     82  ! aerosol : alpha_r=0.001, gas 0.001  (Pruppacher & Klett 1967)
     83  REAL, SAVE :: alpha_r  !  coefficient d'impaction pour la pluie
     84  REAL, SAVE :: alpha_s  !  coefficient d'impaction pour la neige
     85  REAL, SAVE :: R_r      !  mean raindrop radius (m)
     86  REAL, SAVE :: R_s      !  mean snow crystal radius (m)
     87  !$OMP THREADPRIVATE(alpha_r, alpha_s, R_r, R_s)
     88  REAL :: pr, ps, ice, water
    9089  REAL :: conserv
    9190
    92 !!!!!!!!!!!!!!!!!!!! choix lessivage !!!!!!!!!!!!!!!!!!!!!!!!
    93 !!  logical,save  :: inscav_fisrt
    94 !!! $OMP THREADPRIVATE(inscav_first)
    95 
    96 !!!!!!!!!!!!!!!!!!!!!!!!!!!
     91  !!!!!!!!!!!!!!!!!!!! choix lessivage !!!!!!!!!!!!!!!!!!!!!!!!
     92  !!  logical,save  :: inscav_fisrt
     93  !!! $OMP THREADPRIVATE(inscav_first)
     94
     95  !!!!!!!!!!!!!!!!!!!!!!!!!!!
    9796  IF (debut) THEN
    9897
    99 !  inscav_fisrt=.TRUE.
    100 !  CALL getin('inscav_fisrt',inscav_fisrt)
    101 !  IF(inscav_fisrt) THEN
    102 !   PRINT*,'beta from fisrtilp.F90, beta = (z_cond - z_oliq)/z_cond, inscav_fisrt=',inscav_fisrt
    103 !  else
    104 !   PRINT*,'beta from Reddy and Bocuher 2004 (original version), inscav_fisrt=',inscav_fisrt
    105 !  ENDIF
    106 
    107       alpha_r=0.001        !  coefficient d'impaction pour la pluie
    108       alpha_s=0.01         !  coefficient d'impaction pour la neige 
    109       R_r=0.001            !  mean raindrop radius (m)
    110       R_s=0.001            !  mean snow crystal radius (m)
    111       frac_fine_scav=0.7
    112       frac_coar_scav=0.7
    113 !     frac_aer=0.5 ~ droplet size shrinks by evap
    114       frac_aer=0.5
    115 
    116 !JE to speed up, commented 20140219
    117 
    118 !      OPEN(99,file='lsc_scav_param.data',status='old', &
    119 !                form='formatted',err=9999)
    120 !      READ(99,*,end=9998)  alpha_r
    121 !      READ(99,*,end=9998)  alpha_s
    122 !      READ(99,*,end=9998)  R_r
    123 !      READ(99,*,end=9998)  R_s
    124 !      READ(99,*,end=9998)  frac_fine_scav
    125 !      READ(99,*,end=9998)  frac_coar_scav
    126 !      READ(99,*,end=9998)  frac_aer
    127 !9998  Continue
    128 !      CLOSE(99)
    129 !9999  Continue
    130 
    131 !   PRINT*,'alpha_r',alpha_r
    132 !   PRINT*,'alpha_s',alpha_s
    133 !   PRINT*,'R_r',R_r
    134 !   PRINT*,'R_s',R_s
    135 !   PRINT*,'frac_fine_scav',frac_fine_scav
    136 !   PRINT*,'frac_coar_scav',frac_coar_scav
    137 !   PRINT*,'frac_aer ev',frac_aer
    138 
    139 ! JE endcomment
     98    !  inscav_fisrt=.TRUE.
     99    !  CALL getin('inscav_fisrt',inscav_fisrt)
     100    !  IF(inscav_fisrt) THEN
     101    !   PRINT*,'beta from fisrtilp.F90, beta = (z_cond - z_oliq)/z_cond, inscav_fisrt=',inscav_fisrt
     102    !  else
     103    !   PRINT*,'beta from Reddy and Bocuher 2004 (original version), inscav_fisrt=',inscav_fisrt
     104    !  ENDIF
     105
     106    alpha_r = 0.001        !  coefficient d'impaction pour la pluie
     107    alpha_s = 0.01         !  coefficient d'impaction pour la neige
     108    R_r = 0.001            !  mean raindrop radius (m)
     109    R_s = 0.001            !  mean snow crystal radius (m)
     110    frac_fine_scav = 0.7
     111    frac_coar_scav = 0.7
     112    !     frac_aer=0.5 ~ droplet size shrinks by evap
     113    frac_aer = 0.5
     114
     115    !JE to speed up, commented 20140219
     116
     117    !      OPEN(99,file='lsc_scav_param.data',status='old', &
     118    !                form='formatted',err=9999)
     119    !      READ(99,*,end=9998)  alpha_r
     120    !      READ(99,*,end=9998)  alpha_s
     121    !      READ(99,*,end=9998)  R_r
     122    !      READ(99,*,end=9998)  R_s
     123    !      READ(99,*,end=9998)  frac_fine_scav
     124    !      READ(99,*,end=9998)  frac_coar_scav
     125    !      READ(99,*,end=9998)  frac_aer
     126    !9998  Continue
     127    !      CLOSE(99)
     128    !9999  Continue
     129
     130    !   PRINT*,'alpha_r',alpha_r
     131    !   PRINT*,'alpha_s',alpha_s
     132    !   PRINT*,'R_r',R_r
     133    !   PRINT*,'R_s',R_s
     134    !   PRINT*,'frac_fine_scav',frac_fine_scav
     135    !   PRINT*,'frac_coar_scav',frac_coar_scav
     136    !   PRINT*,'frac_aer ev',frac_aer
     137
     138    ! JE endcomment
    140139
    141140  ENDIF !(debut)
    142 !!!!!!!!!!!!!!!!!!!!!!!!!!!
    143 
    144 ! initialization
    145   dxin=0.
    146   dxev=0.
    147   beta_ev=0.
    148 
    149   DO i=1,klon
    150    his_dh(i)=0.
    151   ENDDO
    152 
    153   DO k=1,klev
    154    DO i=1, klon
    155     dxbc(i,k)=0.
    156     beta_v1(i,k)=0.
    157     deltaP(i,k)=0.
    158    ENDDO
    159   ENDDO
    160 
    161   DO k=1,klev
    162     DO i=1, klon
    163      d_tr_insc(i,k,it)=0.
    164      d_tr_bcscav(i,k,it)=0.
    165      d_tr_evap(i,k,it)=0.
    166     ENDDO
    167   ENDDO
    168 
    169 !  pressure and size of the layer
    170   DO k=klev-1, 1, -1
    171    DO i=1, klon
    172      zrho(i,k)=pplay(i,k)/t(i,k)/RD   
    173      zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
    174      zmass(i,k)=(paprs(i,k)-paprs(i,k+1))/RG
    175    ENDDO
    176   ENDDO
    177 
    178     IF (it>1) THEN                               !  aerosol
    179       frac_ev=frac_aer
    180     ELSE                                                !  gas
    181       frac_ev=frac_gas
    182     ENDIF
    183 
    184     IF(it>1) then  ! aerosol
    185      DO k=1, klev
    186       DO i=1, klon
    187        scav(i,k)=frac_fine_scav
    188       ENDDO
    189      ENDDO
    190     ELSE                  ! gas
    191      DO k=1, klev
    192       DO i=1, klon
    193        henry_t=henry*exp(-kk*(1./298.-1./t(i,k)))    !  mol/l/atm
    194        f_a=henry_t/101.325*R*t(i,k)*oliq*zrho(i,k)/rho_water
    195        scav(i,k)=f_a/(1.+f_a)
    196       ENDDO
    197      ENDDO
    198     ENDIF
    199 
    200    DO k=klev-1, 1, -1
    201     DO i=1, klon
    202 !  incloud scavenging
    203 !   IF(inscav_fisrt) THEN
    204    IF (iflag_lscav == 4) THEN
    205       beta=beta_fisrt(i,k)*rneb(i,k)
    206    else
    207       beta=flxr(i,k)-flxr(i,k+1)+flxs(i,k)-flxs(i,k+1)
    208 !      beta=beta/zdz(i,k)/oliq/zrho(i,k)
    209       beta=beta/zmass(i,k)/oliq
    210       beta=MAX(0.,beta)
    211    endif ! (iflag_lscav .EQ. 4)
    212    beta_v1(i,k)=beta    !! for output
    213 
    214       dxin=tr_seri(i,k,it)*(exp(-scav(i,k)*beta*pdtime)-1.)
    215 !      his_dh(i)=his_dh(i)-dxin*zrho(i,k)*zdz(i,k)/pdtime !  kg/m2/s
    216       his_dh(i)=his_dh(i)-dxin*zmass(i,k)/pdtime !  kg/m2/s
    217       d_tr_insc(i,k,it)=dxin
    218 
    219 !  below-cloud impaction
    220     IF(it==1) THEN
    221       d_tr_bcscav(i,k,it)=0.
    222     ELSE
    223      pr=0.5*(flxr(i,k)+flxr(i,k+1))
    224      ps=0.5*(flxs(i,k)+flxs(i,k+1))
    225      water=pr*alpha_r/R_r/rho_water
    226      ice=ps*alpha_s/R_s/rho_ice
    227      dxbc(i,k)=-3./4.*tr_seri(i,k,it)*pdtime*(water+ice)
    228 !   add tracers from below cloud scav in his_dh
    229      his_dh(i)=his_dh(i)-dxbc(i,k)*zmass(i,k)/pdtime !  kg/m2/s
    230      d_tr_bcscav(i,k,it)=dxbc(i,k)
    231     ENDIF
    232 
    233 !  reevaporation
    234       deltaP(i,k)=flxr(i,k+1)+flxs(i,k+1)-flxr(i,k)-flxs(i,k)
    235       deltaP(i,k)=max(deltaP(i,k),0.)
    236 
    237       IF(flxr(i,k+1)+flxs(i,k+1)>1.e-16) THEN
    238        beta_ev(i,k)=deltaP(i,k)/(flxr(i,k+1)+flxs(i,k+1))
     141  !!!!!!!!!!!!!!!!!!!!!!!!!!!
     142
     143  ! initialization
     144  dxin = 0.
     145  dxev = 0.
     146  beta_ev = 0.
     147
     148  DO i = 1, klon
     149    his_dh(i) = 0.
     150  ENDDO
     151
     152  DO k = 1, klev
     153    DO i = 1, klon
     154      dxbc(i, k) = 0.
     155      beta_v1(i, k) = 0.
     156      deltaP(i, k) = 0.
     157    ENDDO
     158  ENDDO
     159
     160  DO k = 1, klev
     161    DO i = 1, klon
     162      d_tr_insc(i, k, it) = 0.
     163      d_tr_bcscav(i, k, it) = 0.
     164      d_tr_evap(i, k, it) = 0.
     165    ENDDO
     166  ENDDO
     167
     168  !  pressure and size of the layer
     169  DO k = klev - 1, 1, -1
     170    DO i = 1, klon
     171      zrho(i, k) = pplay(i, k) / t(i, k) / RD
     172      zdz(i, k) = (paprs(i, k) - paprs(i, k + 1)) / zrho(i, k) / RG
     173      zmass(i, k) = (paprs(i, k) - paprs(i, k + 1)) / RG
     174    ENDDO
     175  ENDDO
     176
     177  IF (it>1) THEN                               !  aerosol
     178    frac_ev = frac_aer
     179  ELSE                                                !  gas
     180    frac_ev = frac_gas
     181  ENDIF
     182
     183  IF(it>1) then  ! aerosol
     184    DO k = 1, klev
     185      DO i = 1, klon
     186        scav(i, k) = frac_fine_scav
     187      ENDDO
     188    ENDDO
     189  ELSE                  ! gas
     190    DO k = 1, klev
     191      DO i = 1, klon
     192        henry_t = henry * exp(-kk * (1. / 298. - 1. / t(i, k)))    !  mol/l/atm
     193        f_a = henry_t / 101.325 * R * t(i, k) * oliq * zrho(i, k) / rho_water
     194        scav(i, k) = f_a / (1. + f_a)
     195      ENDDO
     196    ENDDO
     197  ENDIF
     198
     199  DO k = klev - 1, 1, -1
     200    DO i = 1, klon
     201      !  incloud scavenging
     202      !   IF(inscav_fisrt) THEN
     203      IF (iflag_lscav == 4) THEN
     204        beta = beta_fisrt(i, k) * rneb(i, k)
    239205      else
    240        beta_ev(i,k)=0.
     206        beta = flxr(i, k) - flxr(i, k + 1) + flxs(i, k) - flxs(i, k + 1)
     207        !      beta=beta/zdz(i,k)/oliq/zrho(i,k)
     208        beta = beta / zmass(i, k) / oliq
     209        beta = MAX(0., beta)
     210      endif ! (iflag_lscav .EQ. 4)
     211      beta_v1(i, k) = beta    !! for output
     212
     213      dxin = tr_seri(i, k, it) * (exp(-scav(i, k) * beta * pdtime) - 1.)
     214      !      his_dh(i)=his_dh(i)-dxin*zrho(i,k)*zdz(i,k)/pdtime !  kg/m2/s
     215      his_dh(i) = his_dh(i) - dxin * zmass(i, k) / pdtime !  kg/m2/s
     216      d_tr_insc(i, k, it) = dxin
     217
     218      !  below-cloud impaction
     219      IF(it==1) THEN
     220        d_tr_bcscav(i, k, it) = 0.
     221      ELSE
     222        pr = 0.5 * (flxr(i, k) + flxr(i, k + 1))
     223        ps = 0.5 * (flxs(i, k) + flxs(i, k + 1))
     224        water = pr * alpha_r / R_r / rho_water
     225        ice = ps * alpha_s / R_s / rho_ice
     226        dxbc(i, k) = -3. / 4. * tr_seri(i, k, it) * pdtime * (water + ice)
     227        !   add tracers from below cloud scav in his_dh
     228        his_dh(i) = his_dh(i) - dxbc(i, k) * zmass(i, k) / pdtime !  kg/m2/s
     229        d_tr_bcscav(i, k, it) = dxbc(i, k)
     230      ENDIF
     231
     232      !  reevaporation
     233      deltaP(i, k) = flxr(i, k + 1) + flxs(i, k + 1) - flxr(i, k) - flxs(i, k)
     234      deltaP(i, k) = max(deltaP(i, k), 0.)
     235
     236      IF(flxr(i, k + 1) + flxs(i, k + 1)>1.e-16) THEN
     237        beta_ev(i, k) = deltaP(i, k) / (flxr(i, k + 1) + flxs(i, k + 1))
     238      else
     239        beta_ev(i, k) = 0.
    241240      endif
    242241
    243       beta_ev(i,k)=max(min(1.,beta_ev(i,k)),0.)
    244 
    245 !jyg
    246      
    247       IF(abs(1-(1-frac_ev)*beta_ev(i,k))>1.e-16) THEN
    248 ! remove tracers from precipitation owing to release by evaporation in his_dh
    249 !      dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
    250       dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zmass(i,k)) &
    251                                       /(1 -(1-frac_ev)*beta_ev(i,k))
    252        his_dh(i)=his_dh(i)*(1 - frac_ev*beta_ev(i,k) / (1 -(1-frac_ev)*beta_ev(i,k)))
     242      beta_ev(i, k) = max(min(1., beta_ev(i, k)), 0.)
     243
     244      !jyg
     245
     246      IF(abs(1 - (1 - frac_ev) * beta_ev(i, k))>1.e-16) THEN
     247        ! remove tracers from precipitation owing to release by evaporation in his_dh
     248        !      dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
     249        dxev = frac_ev * beta_ev(i, k) * his_dh(i) * pdtime / (zmass(i, k)) &
     250                / (1 - (1 - frac_ev) * beta_ev(i, k))
     251        his_dh(i) = his_dh(i) * (1 - frac_ev * beta_ev(i, k) / (1 - (1 - frac_ev) * beta_ev(i, k)))
    253252      else
    254 !       dxev=his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k))
    255        dxev=his_dh(i) *pdtime/(zmass(i,k))
    256        his_dh(i)=0.
     253        !       dxev=his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k))
     254        dxev = his_dh(i) * pdtime / (zmass(i, k))
     255        his_dh(i) = 0.
    257256      endif
    258 !      PRINT*,  k, 'beta_ev',beta_ev
    259 ! remove tracers from precipitation owing to release by evaporation in his_dh
    260 !!      dxev=frac_ev*deltaP(i,k)*pdtime * his_dh(i) /(zrho(i,k)*zdz(i,k))
    261 !rplmd
    262 !!      dxev=frac_ev*deltaP(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
    263 !!                                      /max(flxr(i,k)+flxs(i,k),1.e-16)
    264 
    265 
    266       d_tr_evap(i,k,it)=dxev
    267 !!     tendency is further added in phytrac x = x + dx
     257      !      PRINT*,  k, 'beta_ev',beta_ev
     258      ! remove tracers from precipitation owing to release by evaporation in his_dh
     259      !!      dxev=frac_ev*deltaP(i,k)*pdtime * his_dh(i) /(zrho(i,k)*zdz(i,k))
     260      !rplmd
     261      !!      dxev=frac_ev*deltaP(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
     262      !!                                      /max(flxr(i,k)+flxs(i,k),1.e-16)
     263
     264      d_tr_evap(i, k, it) = dxev
     265      !!     tendency is further added in phytrac x = x + dx
    268266    ENDDO !!  do i
    269    ENDDO  !! do k
    270 
    271 !jyg (20130114)
    272    DO i = 1,klon
    273      qPrls(i,it) = his_dh(i)/max(flxr(i,1)+flxs(i,1),1.e-16)
    274    ENDDO
    275 !jyg end
    276 
    277 
    278 ! test de conservation
    279       conserv=0.
    280 !      DO k= klev,1,-1
    281 !        DO i=1, klon
    282 !         conserv=conserv+d_tr_insc(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG &
    283 !                +d_tr_bcscav(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG  &
    284 !                +d_tr_evap(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG
    285 !      IF(it.EQ.3) WRITE(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)'),&
    286 !      k,'lsc conserv ',conserv,'insc',d_tr_insc(i,k,it),'bc',d_tr_bcscav(i,k,it),'ev',d_tr_evap(i,k,it)
    287 !       ENDDO
    288 !     ENDDO
     267  ENDDO  !! do k
     268
     269  !jyg (20130114)
     270  DO i = 1, klon
     271    qPrls(i, it) = his_dh(i) / max(flxr(i, 1) + flxs(i, 1), 1.e-16)
     272  ENDDO
     273  !jyg end
     274
     275
     276  ! test de conservation
     277  conserv = 0.
     278  !      DO k= klev,1,-1
     279  !        DO i=1, klon
     280  !         conserv=conserv+d_tr_insc(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG &
     281  !                +d_tr_bcscav(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG  &
     282  !                +d_tr_evap(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG
     283  !      IF(it.EQ.3) WRITE(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)'),&
     284  !      k,'lsc conserv ',conserv,'insc',d_tr_insc(i,k,it),'bc',d_tr_bcscav(i,k,it),'ev',d_tr_evap(i,k,it)
     285  !       ENDDO
     286  !     ENDDO
    289287
    290288END SUBROUTINE lsc_scav_orig
Note: See TracChangeset for help on using the changeset viewer.