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

Last change on this file since 2312 was 2312, checked in by lrossi, 5 years ago

MARS GCM:
Implementation of HDO in the physics
New tracers hdo_vap and hdo_ice are added. Cf. traceur.def.isotopes in deftank for exemple of traceur.def.
All HDO related computations are activated by flag hdo=.true. in callphys.def. (see callphys.def.hdo in deftank for an example)
Additional option hdofrac (true by default) allows for turning on/off fractionation (for tests).
For now, only functional with simplified cloud scheme (so microphys=.false. and activice=.false. also recommended).
Initialisation of start files can be done with option 'inihdo' in newstart.
LR

File size: 5.2 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,zq,pqsurf,
9     &                      old_h2o_vap,pdqsdif,h2oflux,
10     &                      hdoflux)
11
12      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice,
13     &                      igcm_hdo_vap, igcm_hdo_ice
14
15      use surfdat_h, only: watercaptag
16      use geometry_mod, only: longitude_deg,latitude_deg
17
18      implicit none
19c------------------------------------------------------------------
20c               Routine to compute the fluxes between air and surface
21c               for HDO, based of the fluxes for H2O
22c           L. Rossi.; M. Vals 2019
23c------------------------------------------------------------------
24      include "callkeys.h"
25
26c------------------------------------------------------------------
27c     Arguments:
28c     ---------
29c     Inputs:
30      INTEGER, INTENT(IN) :: ngrid,nlay
31      INTEGER, INTENT(IN) :: nq                 ! nombre de traceurs
32
33      REAL, INTENT(IN) :: ptimestep             ! pas de temps physique (s)
34      REAL, INTENT(IN) :: zt(ngrid,nlay)       ! local value of temperature
35      REAL, INTENT(IN) :: zq(ngrid,nlay,nq)    ! local value of tracers
36      REAL, INTENT(IN) :: pqsurf(ngrid,nq)
37      REAL, INTENT(IN) :: old_h2o_vap(ngrid)     ! traceur d'eau avant
38                                           !traitement de l'eau (kg/kg)
39      REAL, INTENT(INOUT) :: pdqsdif(ngrid,nq)    ! tendance towards surface
40                                 !   (kg/kg.s-1)
41      REAL, INTENT(IN) :: h2oflux(ngrid)       ! value of vapour flux of H2O
42                                      ! same sign as pdqsdif
43
44c     Output:
45      REAL, INTENT(OUT) :: hdoflux(ngrid)       ! value of vapour flux of HDO
46
47c------------------------------------------------------------------
48c     Local variables:
49
50      REAL alpha_c(ngrid)  ! fractionation factor
51      REAL extrasublim ! sublimation in excess of surface ice
52      REAL tmpratio(ngrid)   ! D/H ratio in flux to surf
53
54      INTEGER ig,l
55
56      REAL DoH_vap(ngrid)
57
58c-----------------------------------------------------------------------
59c    Calculation of the fluxes for HDO
60       
61        alpha_c(1:ngrid)=0.
62
63        DO ig=1,ngrid
64
65            !! IF Sublimation
66            if (h2oflux(ig).le.0.) then
67
68               if (pqsurf(ig,igcm_h2o_ice).gt.1e-16) then
69                pdqsdif(ig,igcm_hdo_ice) =
70     &            pdqsdif(ig,igcm_h2o_ice)*
71     &             (pqsurf(ig,igcm_hdo_ice)/
72     &             pqsurf(ig,igcm_h2o_ice) )
73               else
74                pdqsdif(ig,igcm_hdo_ice) = 0.
75               endif
76
77                pdqsdif(ig,igcm_hdo_ice)=
78     &             max(pdqsdif(ig,igcm_hdo_ice),
79     &            -pqsurf(ig,igcm_hdo_ice)/ptimestep)
80
81                hdoflux(ig) = pdqsdif(ig,igcm_hdo_ice)
82
83             if(watercaptag(ig)) then
84
85              !if we sublimate more than qsurf
86              if ((-h2oflux(ig)*ptimestep)
87     &           .gt.pqsurf(ig,igcm_h2o_ice)) then
88
89C               This is how much we sublimate in excess of
90C               pqsurf for H2O                       
91                  extrasublim=
92     &              ( pqsurf(ig,igcm_h2o_ice) /ptimestep
93     &                   +h2oflux(ig))
94c               extrasublim has same sign conv. as pdqsdif
95 
96C               hdoflux(ig) is the flux of HDO from atm. to surf.
97c               The D/H of the old ice is supposed to be 5 SMOW
98c               We need D/H of the flux to be 5, so we need
99c               extrasublim* 5 * 2 * 155.76e-6 (=1 SMOW)
100                    hdoflux(ig)= hdoflux(ig)
101     &                       +(extrasublim*(2.*155.76e-6)*5.)
102!    &                       +extrasublim
103!                 hdoflux(ig) = h2oflux(ig) !test
104              endif
105                endif ! watercap
106
107            else ! condensation
108
109               if (hdofrac) then !do we use fractionation?
110                alpha_c(ig) = exp(16288./zt(ig,1)**2.-9.34e-2)
111c               alpha_c = exp(13525./zt(ig,l)**2.-5.59e-2) !Lamb
112               else
113                alpha_c(ig) = 1.
114               endif
115
116               if (old_h2o_vap(ig).gt.1.e-16) then
117                         pdqsdif(ig,igcm_hdo_ice)=
118     &                      alpha_c(ig)*pdqsdif(ig,igcm_h2o_ice)*
119     &                      (zq(ig,1,igcm_hdo_vap)/
120     &                          old_h2o_vap(ig))
121                else
122                 pdqsdif(ig,igcm_hdo_ice)= 0.
123                endif
124
125               if (hdofrac) then !do we use fractionation?
126                pdqsdif(ig,igcm_hdo_ice)=
127     &      min(  pdqsdif(ig,igcm_hdo_ice),
128     &           (zq(ig,1,igcm_hdo_vap)/ptimestep) )
129               endif
130
131                hdoflux(ig)=pdqsdif(ig,igcm_hdo_ice)
132
133            endif !sublimation
134
135        ENDDO ! of DO ig=1,ngrid
136
137c           CALL WRITEDIAGFI(ngrid,'extrasublim',
138c    &                       'extrasublimation',
139c    &                       ' ',2,tmpratio)
140c           CALL WRITEDIAGFI(ngrid,'alpha_c_s',
141c    &                       'alpha_c_s',
142c    &                       ' ',2,alpha_c)
143
144       return
145      end subroutine hdo_surfex
146c------------------------------------------------------------------
147
148      end module hdo_surfex_mod
Note: See TracBrowser for help on using the repository browser.