- Timestamp:
- Jul 23, 2024, 5:57:06 PM (2 months ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/inscav_spl.f90
r5103 r5104 1 SUBROUTINE inscav_spl(pdtime,it,masse,henry,kk,qliq, 2 . flxr,flxs,zrho,zdz,t,x,3 . his_dh)4 5 IMPLICIT NONE6 c=====================================================================7 cObjet : depot humide de traceurs8 cDate : mars 19989 c Auteur: O. Boucher (LOA) 10 c=====================================================================11 c 12 13 14 15 INCLUDE "YOECUMF.h"16 c 17 INTEGERit18 REALpdtime ! pas de temps (s)19 REALmasse ! molar mass (except for BC/OM/IF/DUST=Nav)20 REALhenry ! constante de Henry en mol/l/atm21 REALkk ! coefficient de dependence en T (K)22 REALqliq ! contenu en eau liquide dans le nuage (kg/kg)23 !REAL flxr(klon,klev+1) ! flux precipitant de pluie24 !REAL flxs(klon,klev+1) ! flux precipitant de neige25 REAL flxr(klon,klev) ! flux precipitant de pluie ! Titane26 REAL flxs(klon,klev) ! flux precipitant de neige ! Titane27 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 traceur33 REALhis_dh(klon) ! tendance de traceur integre verticalement34 c 35 c--variables locales 36 INTEGERi, k37 c 38 REALdx ! tendance de traceur39 REALf_a !--rapport de la phase aqueuse a la phase gazeuse40 REALbeta !--taux de conversion de l'eau en pluie41 REALhenry_t !--constante de Henry a T t (mol/l/atm)42 REAL scav(klon,klev) !--fraction aqueuse du constituant43 REAL K1, K2, ph, frac44 REAL frac_gas, frac_aer !-cste pour la reevaporation45 PARAMETER (ph=5., frac_gas=1.0, frac_aer=0.5)46 c---cste de dissolution pour le depot humide47 REAL frac_fine_scav,frac_coar_scav48 c---added by nhl49 REALaux_cte1 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 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 50 51 PARAMETER (frac_fine_scav=0.7)52 PARAMETER (frac_coar_scav=0.7)51 PARAMETER (frac_fine_scav = 0.7) 52 PARAMETER (frac_coar_scav = 0.7) 53 53 54 c--101.325 m3/l x Pa/atm 55 c--R Pa.m3/mol/K 56 c 57 c------------------------------------------ 58 c 59 cnhl 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 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) 64 127 ENDIF 65 c 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 c 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 c--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 c--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 c 139 RETURN 140 END 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 140 END SUBROUTINE inscav_spl
Note: See TracChangeset
for help on using the changeset viewer.
