Context Navigation

hdo_surfex_mod.F @ 3799

Last change on this file since 3799 was 3726, checked in by emillour, 2 months ago
Mars PCM: Turn "callkeys.h" into module "callkeys_mod.F90" EM
File size: 6.7 KB

Line
1	MODULE hdo_surfex_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine hdo_surfex(ngrid,nlay,nq,ptimestep,
8	& zt,pplay,zq,pqsurf,
9	& old_h2o_vap,qsat,pdqsdif,dwatercap_dif,
10	& hdoflux)
11
12	use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice,
13	& igcm_hdo_vap, igcm_hdo_ice,
14	& qperemin
15	use surfdat_h, only: watercaptag
16	use geometry_mod, only: longitude_deg,latitude_deg
17	use comcstfi_h, only: pi
18	use microphys_h, only: nav, kbz, mh2o, mco2, mhdo
19	use microphys_h, only: molco2, molh2o, molhdo
20	use write_output_mod, only: write_output
21	use callkeys_mod, only: hdofrac
22
23	implicit none
24	c------------------------------------------------------------------
25	c Routine to compute the fluxes between air and surface
26	c for HDO, based of the fluxes for H2O
27	c L. Rossi.; M. Vals 2019
28	c------------------------------------------------------------------
29	c Arguments:
30	c ---------
31	c Inputs:
32	INTEGER, INTENT(IN) :: ngrid,nlay
33	INTEGER, INTENT(IN) :: nq ! nombre de traceurs
34
35	REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s)
36	REAL, INTENT(IN) :: zt(ngrid,nlay) ! local value of temperature
37	REAL, INTENT(IN) :: pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
38	REAL, INTENT(IN) :: zq(ngrid,nlay,nq) ! local value of tracers
39	REAL, INTENT(IN) :: pqsurf(ngrid,nq)
40	REAL, INTENT(IN) :: old_h2o_vap(ngrid) ! traceur d'eau avant
41	!traitement de l'eau (kg/kg)
42	REAL, INTENT(IN) :: qsat(ngrid) ! saturation mixing ratio
43	REAL, INTENT(IN) :: dwatercap_dif(ngrid) ! trend related to permanent ice
44	REAL, INTENT(INOUT) :: pdqsdif(ngrid,nq) ! tendance towards surface
45	! (kg/kg.s-1)
46
47	c Output:
48	REAL, INTENT(OUT) :: hdoflux(ngrid) ! value of vapour flux of HDO
49
50	c------------------------------------------------------------------
51	c Local variables:
52
53	REAL alpha(ngrid) ! equilibrium fractionation factor
54	REAL alpha_c(ngrid) ! real fractionation factor
55	REAL extrasublim ! sublimation in excess of surface ice
56	REAL tmpratio(ngrid) ! D/H ratio in flux to surf
57	REAL h2oflux(ngrid) ! value of vapour flux of H2O
58	! same sign as pdqsdif
59	REAL*8 satu(ngrid) ! Water vapor saturation ratio over ice
60	REAL zt1(ngrid),pplay1(ngrid)
61	REAL Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
62	INTEGER ig,l
63
64	REAL DoH_vap(ngrid)
65
66	c-----------------------------------------------------------------------
67	c Calculation of the fluxes for HDO
68	!! Calculation of the saturation ratio in the layer above the surface
69	satu(1:ngrid)=old_h2o_vap(1:ngrid) / qsat(1:ngrid)
70	!! Initialisation of the fractionation coefficient
71	alpha(1:ngrid)=1.
72	alpha_c(1:ngrid)=1.
73
74	DO ig=1,ngrid
75
76	h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice) +
77	& dwatercap_dif(ig)
78
79	!! IF Sublimation
80	if (h2oflux(ig).le.0.) then
81
82	if (pqsurf(ig,igcm_h2o_ice).gt.qperemin) then
83	pdqsdif(ig,igcm_hdo_ice) =
84	& pdqsdif(ig,igcm_h2o_ice)*
85	& (pqsurf(ig,igcm_hdo_ice)/
86	& pqsurf(ig,igcm_h2o_ice) )
87	else
88	pdqsdif(ig,igcm_hdo_ice) = 0.
89	endif
90
91	pdqsdif(ig,igcm_hdo_ice)=
92	& max(pdqsdif(ig,igcm_hdo_ice),
93	& -pqsurf(ig,igcm_hdo_ice)/ptimestep)
94
95	hdoflux(ig) = pdqsdif(ig,igcm_hdo_ice)
96
97	if(watercaptag(ig)) then
98
99	!if we sublimate more than qsurf
100	if ((-h2oflux(ig)*ptimestep)
101	& .gt.pqsurf(ig,igcm_h2o_ice)) then
102
103	C dwatercap_dif is how much we sublimate in excess of
104	C pqsurf for H2O
105	C hdoflux(ig) is the flux of HDO from atm. to surf.
106	c The D/H of the old ice is supposed to be 5 SMOW
107	c We need D/H of the flux to be 5, so we need
108	c dwatercap_dif* 5 * 2 * 155.76e-6 (=1 SMOW)
109	hdoflux(ig)= hdoflux(ig)
110	& +(dwatercap_dif(ig)(2.155.76e-6)*5.)
111	endif
112	endif ! watercap
113
114	else ! condensation
115
116	if (hdofrac) then !do we use fractionation?
117	!! Calculation of the H2O vapor diffusion coefficient
118	Dv = 1./3. * sqrt( 8kbzzt(ig,1)/(pi*mh2o/nav) )
119	& * kbz * zt(ig,1) /
120	& ( pi * pplay(ig,1) * (molco2+molh2o)*(molco2+molh2o)
121	& * sqrt(1.+mh2o/mco2) )
122	!! Calculation of the HDO vapor diffusion coefficient
123	Dv_hdo = 1./3. * sqrt( 8kbzzt(ig,1)/(pi*mhdo/nav) )
124	& * kbz * zt(ig,1) /
125	& ( pi * pplay(ig,1) * (molco2+molhdo)*(molco2+molhdo)
126	& * sqrt(1.+mhdo/mco2) )
127	!! Calculation of the "equilibrium" fractionation coefficient
128	c alpha(ig) = exp(16288./zt(ig,1)**2.-9.34e-2)
129	alpha(ig) = exp(13525./zt(ig,1)**2.-5.59e-2) !Lamb
130	!! Calculation of the 'real' fractionnation coefficient (effect of kinetics, see Jouzel and Merlivat, 1984)
131	alpha_c(ig) = (alpha(ig)*satu(ig))/
132	& ( (alpha(ig)(Dv/Dv_hdo)(satu(ig)-1.)) + 1.)
133	else
134	alpha_c(ig) = 1.
135	endif
136
137	if (old_h2o_vap(ig).gt.qperemin) then
138	pdqsdif(ig,igcm_hdo_ice)=
139	& alpha_c(ig)pdqsdif(ig,igcm_h2o_ice)
140	& (zq(ig,1,igcm_hdo_vap)/
141	& old_h2o_vap(ig))
142	else
143	pdqsdif(ig,igcm_hdo_ice)= 0.
144	endif
145
146	if (hdofrac) then !do we use fractionation?
147	pdqsdif(ig,igcm_hdo_ice)=
148	& min( pdqsdif(ig,igcm_hdo_ice),
149	& (zq(ig,1,igcm_hdo_vap)/ptimestep) )
150	endif
151
152	hdoflux(ig)=pdqsdif(ig,igcm_hdo_ice)
153
154	endif !sublimation
155
156	ENDDO ! of DO ig=1,ngrid
157
158	c CALL write_output('extrasublim',
159	c & 'extrasublimation',
160	c & ' ',tmpratio)
161	c CALL write_output('alpha_c_s',
162	c & 'alpha_c_s',
163	c & ' ',alpha_c)
164
165	end subroutine hdo_surfex
166	c------------------------------------------------------------------
167
168	end module hdo_surfex_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format