source: LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/inscav_spl.f90 @ 5104

Last change on this file since 5104 was 5104, checked in by abarral, 2 months ago

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F in DUST to *.f90
File size: 4.9 KB
Line 
1SUBROUTINE inscav_spl(pdtime, it, masse, henry, kk, qliq, &
2        flxr, flxs, zrho, zdz, t, x, &
3        his_dh)
4  USE dimphy
5  IMPLICIT NONE
6  !=====================================================================
7  ! Objet : depot humide de traceurs
8  ! Date : mars 1998
9  ! Auteur: O. Boucher (LOA)
10  !=====================================================================
11  !
12  INCLUDE "dimensions.h"
13  INCLUDE "chem.h"
14  INCLUDE "YOMCST.h"
15  INCLUDE "YOECUMF.h"
16  !
17  INTEGER :: it
18  REAL :: pdtime              ! pas de temps (s)
19  REAL :: masse               ! molar mass (except for BC/OM/IF/DUST=Nav)
20  REAL :: henry               ! constante de Henry en mol/l/atm
21  REAL :: kk                  ! coefficient de dependence en T (K)
22  REAL :: qliq                ! contenu en eau liquide dans le nuage (kg/kg)
23  ! REAL flxr(klon,klev+1)   ! flux precipitant de pluie
24  ! REAL flxs(klon,klev+1)   ! flux precipitant de neige
25  REAL :: flxr(klon, klev)   ! flux precipitant de pluie   ! Titane
26  REAL :: flxs(klon, klev)   ! flux precipitant de neige   ! Titane
27  REAL :: flxr_aux(klon, klev + 1)
28  REAL :: flxs_aux(klon, klev + 1)
29  REAL :: zrho(klon, klev)
30  REAL :: zdz(klon, klev)
31  REAL :: t(klon, klev)
32  REAL :: x(klon, klev)        ! q de traceur
33  REAL :: his_dh(klon)        ! tendance de traceur integre verticalement
34  !
35  !--variables locales
36  INTEGER :: i, k
37  !
38  REAL :: dx      ! tendance de traceur
39  REAL :: f_a     !--rapport de la phase aqueuse a la phase gazeuse
40  REAL :: beta    !--taux de conversion de l'eau en pluie
41  REAL :: henry_t !--constante de Henry a T t  (mol/l/atm)
42  REAL :: scav(klon, klev)    !--fraction aqueuse du constituant
43  REAL :: K1, K2, ph, frac
44  REAL :: frac_gas, frac_aer !-cste pour la reevaporation
45  PARAMETER (ph = 5., frac_gas = 1.0, frac_aer = 0.5)
46  !---cste de dissolution pour le depot humide
47  REAL :: frac_fine_scav, frac_coar_scav
48  !---added by nhl
49  REAL :: aux_cte
50
51  PARAMETER (frac_fine_scav = 0.7)
52  PARAMETER (frac_coar_scav = 0.7)
53
54  !--101.325  m3/l x Pa/atm
55  !--R        Pa.m3/mol/K
56  !
57  !------------------------------------------
58  !
59  !nhl      IF (it.EQ.2.OR.it.EQ.3) THEN !--aerosol  ! AS IT WAS FIRST
60  IF (it==2.OR.it==3.OR.it==4) THEN !--aerosol
61    frac = frac_aer
62  ELSE                                                !--gas
63    frac = frac_gas
64  ENDIF
65  !
66  IF (it==1) THEN
67    DO k = 1, klev
68      DO i = 1, klon
69        henry_t = henry * exp(-kk * (1. / 298. - 1. / t(i, k)))    !--mol/l/atm
70        K1 = 1.2e-2 * exp(-2010 * (1 / 298. - 1 / t(i, k)))
71        K2 = 6.6e-8 * exp(-1510 * (1 / 298. - 1 / t(i, k)))
72        henry_t = henry_t * (1 + K1 / 10.**(-ph) + K1 * K2 / (10.**(-ph))**2)
73        f_a = henry_t / 101.325 * R * t(i, k) * qliq * zrho(i, k) / rho_water
74        scav(i, k) = f_a / (1. + f_a)
75      ENDDO
76    ENDDO
77  ELSEIF (it==2) THEN
78    DO k = 1, klev
79      DO i = 1, klon
80        scav(i, k) = frac_fine_scav
81      ENDDO
82    ENDDO
83  ELSEIF (it==3) THEN
84    DO k = 1, klev
85      DO i = 1, klon
86        scav(i, k) = frac_coar_scav
87      ENDDO
88    ENDDO
89  ELSEIF (it==4) THEN
90    DO k = 1, klev
91      DO i = 1, klon
92        scav(i, k) = frac_coar_scav
93      ENDDO
94    ENDDO
95  ELSE
96    PRINT *, 'it non pris en compte'
97    STOP
98  ENDIF
99  !
100  ! NHL
101  ! Auxiliary variables defined to deal with the fact that precipitation
102  ! fluxes are defined on klev levels only.
103  ! NHL
104
105  flxr_aux(:, klev + 1) = 0.0
106  flxs_aux(:, klev + 1) = 0.0
107  flxr_aux(:, 1:klev) = flxr(:, :)
108  flxs_aux(:, 1:klev) = flxs(:, :)
109  DO k = klev, 1, -1
110    DO i = 1, klon
111      !--scavenging
112      beta = flxr_aux(i, k) - flxr_aux(i, k + 1) + flxs_aux(i, k) - flxs_aux(i, k + 1)
113      beta = beta / zdz(i, k) / qliq / zrho(i, k)
114      beta = MAX(0.0, beta)
115      dx = x(i, k) * (exp(-scav(i, k) * beta * pdtime) - 1.)
116      x(i, k) = x(i, k) + dx
117      his_dh(i) = his_dh(i) - dx / RNAVO * &
118              masse * 1.e3 * 1.e6 * zdz(i, k) / pdtime !--mgS/m2/s
119      !--reevaporation
120      beta = flxr_aux(i, k) - flxr_aux(i, k + 1) + flxs_aux(i, k) - flxs_aux(i, k + 1)
121      IF (beta<0.) beta = beta / (flxr_aux(i, k + 1) + flxs_aux(i, k + 1))
122      IF (flxr_aux(i, k) + flxs_aux(i, k)==0) THEN  !--reevaporation totale
123        beta = MIN(MAX(0.0, -beta), 1.0)
124      ELSE                          !--reevaporation non totale pour aerosols
125        ! !print *,'FRAC USED IN INSCAV_SPL'
126        beta = MIN(MAX(0.0, -beta) * frac, 1.0)
127      ENDIF
128      dx = beta * his_dh(i) * RNAVO / masse / 1.e3 / 1.e6 / zdz(i, k) * pdtime !ORIG LINE
129      ! funny line for TL/AD
130      ! AD test works without (x) and for xd = dxd*1.e5 : 2.79051851638 times the 0.
131      ! AD test does not work with the line : 754592404.083 times the 0.
132      ! problem seems to be linked to the largest dx wrt x
133      ! x(i, k) = x(i, k) + dx
134      !  x(i, k) = x(i, k) + dx         ! THIS LINE WAS COMMENTED OUT ORIGINALY !nhl
135      his_dh(i) = (1. - beta) * his_dh(i)
136    ENDDO
137  ENDDO
138  !
139  RETURN
140END SUBROUTINE inscav_spl
Note: See TracBrowser for help on using the repository browser.