source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/inscav_spl.F @ 4934

Last change on this file since 4934 was 2175, checked in by jescribano, 10 years ago

SPLA code included for first time

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