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

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

Mars GCM:
Fixing some errors in vdifc_mod related to variable watercap. This variable was also integrated to the hdo cycle.
Also added watercap output for the 1D model.
LR

File size: 5.1 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,dwatercap_dif,
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(IN) :: dwatercap_dif(ngrid)  ! trend related to permanent ice
40      REAL, INTENT(INOUT) :: pdqsdif(ngrid,nq)    ! tendance towards surface
41                                 !   (kg/kg.s-1)
42
43c     Output:
44      REAL, INTENT(OUT) :: hdoflux(ngrid)       ! value of vapour flux of HDO
45
46c------------------------------------------------------------------
47c     Local variables:
48
49      REAL alpha_c(ngrid)  ! fractionation factor
50      REAL extrasublim ! sublimation in excess of surface ice
51      REAL tmpratio(ngrid)   ! D/H ratio in flux to surf
52      REAL h2oflux(ngrid)       ! value of vapour flux of H2O
53                                      ! same sign as pdqsdif
54
55      INTEGER ig,l
56
57      REAL DoH_vap(ngrid)
58
59c-----------------------------------------------------------------------
60c    Calculation of the fluxes for HDO
61       
62        alpha_c(1:ngrid)=0.
63
64        DO ig=1,ngrid
65             
66            h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice) +
67     &          dwatercap_dif(ig)
68
69            !! IF Sublimation
70            if (h2oflux(ig).le.0.) then
71
72               if (pqsurf(ig,igcm_h2o_ice).gt.1e-16) then
73                pdqsdif(ig,igcm_hdo_ice) =
74     &            pdqsdif(ig,igcm_h2o_ice)*
75     &             (pqsurf(ig,igcm_hdo_ice)/
76     &             pqsurf(ig,igcm_h2o_ice) )
77               else
78                pdqsdif(ig,igcm_hdo_ice) = 0.
79               endif
80
81                pdqsdif(ig,igcm_hdo_ice)=
82     &             max(pdqsdif(ig,igcm_hdo_ice),
83     &            -pqsurf(ig,igcm_hdo_ice)/ptimestep)
84
85                hdoflux(ig) = pdqsdif(ig,igcm_hdo_ice)
86
87             if(watercaptag(ig)) then
88
89              !if we sublimate more than qsurf
90              if ((-h2oflux(ig)*ptimestep)
91     &           .gt.pqsurf(ig,igcm_h2o_ice)) then
92
93C               dwatercap_dif is how much we sublimate in excess of
94C               pqsurf for H2O                       
95C               hdoflux(ig) is the flux of HDO from atm. to surf.
96c               The D/H of the old ice is supposed to be 5 SMOW
97c               We need D/H of the flux to be 5, so we need
98c               dwatercap_dif* 5 * 2 * 155.76e-6 (=1 SMOW)
99                    hdoflux(ig)= hdoflux(ig)
100     &                   +(dwatercap_dif(ig)*(2.*155.76e-6)*5.)
101                endif
102             endif ! watercap
103
104            else ! condensation
105
106               if (hdofrac) then !do we use fractionation?
107                alpha_c(ig) = exp(16288./zt(ig,1)**2.-9.34e-2)
108c               alpha_c = exp(13525./zt(ig,l)**2.-5.59e-2) !Lamb
109               else
110                alpha_c(ig) = 1.
111               endif
112
113               if (old_h2o_vap(ig).gt.1.e-16) then
114                         pdqsdif(ig,igcm_hdo_ice)=
115     &                      alpha_c(ig)*pdqsdif(ig,igcm_h2o_ice)*
116     &                      (zq(ig,1,igcm_hdo_vap)/
117     &                          old_h2o_vap(ig))
118                else
119                 pdqsdif(ig,igcm_hdo_ice)= 0.
120                endif
121
122               if (hdofrac) then !do we use fractionation?
123                pdqsdif(ig,igcm_hdo_ice)=
124     &      min(  pdqsdif(ig,igcm_hdo_ice),
125     &           (zq(ig,1,igcm_hdo_vap)/ptimestep) )
126               endif
127
128                hdoflux(ig)=pdqsdif(ig,igcm_hdo_ice)
129
130            endif !sublimation
131
132        ENDDO ! of DO ig=1,ngrid
133
134c           CALL WRITEDIAGFI(ngrid,'extrasublim',
135c    &                       'extrasublimation',
136c    &                       ' ',2,tmpratio)
137c           CALL WRITEDIAGFI(ngrid,'alpha_c_s',
138c    &                       'alpha_c_s',
139c    &                       ' ',2,alpha_c)
140
141       return
142      end subroutine hdo_surfex
143c------------------------------------------------------------------
144
145      end module hdo_surfex_mod
Note: See TracBrowser for help on using the repository browser.