source: LMDZ6/branches/Ocean_skin/libf/phylmd/Dust/inscav_spl.F @ 3747

Last change on this file since 3747 was 2630, checked in by fhourdin, 8 years ago

Importation du modèle d'aérosols de Boucher, Escribano et al.

File size: 4.7 KB
Line 
1      SUBROUTINE inscav_spl(pdtime,it,masse,henry,kk,qliq,
2     .                   flxr,flxs,zrho,zdz,t,x,
3     .                   his_dh)
4      USE dimphy
5      IMPLICIT NONE
6c=====================================================================
7c Objet : depot humide de traceurs
8c Date : mars 1998
9c Auteur: O. Boucher (LOA)
10c=====================================================================
11c
12#include "dimensions.h"
13#include "chem.h"
14#include "YOMCST.h"
15#include "YOECUMF.h"
16c
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
34c
35c--variables locales     
36      INTEGER i, k
37c
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)
46c---cste de dissolution pour le depot humide
47      REAL frac_fine_scav,frac_coar_scav
48c---added by nhl
49      REAL aux_cte
50
51      PARAMETER (frac_fine_scav=0.7)
52      PARAMETER (frac_coar_scav=0.7)
53
54c--101.325  m3/l x Pa/atm
55c--R        Pa.m3/mol/K
56c
57c------------------------------------------
58c
59cnhl      IF (it.EQ.2.OR.it.EQ.3) THEN !--aerosol  ! AS IT WAS FIRST
60      IF (it.EQ.2.OR.it.EQ.3.OR.it.EQ.4) THEN !--aerosol
61        frac=frac_aer
62      ELSE                                                !--gas
63        frac=frac_gas
64      ENDIF
65c
66      IF (it.EQ.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.EQ.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.EQ.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.EQ.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
99c
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
111c--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
119c--reevaporation
120        beta=flxr_aux(i,k)-flxr_aux(i,k+1)+flxs_aux(i,k)-flxs_aux(i,k+1)
121        IF (beta.LT.0.) beta=beta/(flxr_aux(i,k+1)+flxs_aux(i,k+1))
122        IF (flxr_aux(i,k)+flxs_aux(i,k).EQ.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
138c
139      RETURN
140      END
Note: See TracBrowser for help on using the repository browser.