source: trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F @ 2800

Last change on this file since 2800 was 2593, checked in by mvals, 3 years ago

Mars GCM:
hdo_surfex_mod.F, vdifc_mod.F: addings to account for the effect of kinetics in the fractionation by condensation of HDO at the surface
MV

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