Context Navigation

hdo_surfex_mod.F @ 3400

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

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

Download in other formats:

Original Format