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

Last change on this file since 2390 was 2378, checked in by lrossi, 4 years ago

MARS GCM

HDO
Correction of an error in newstart for inihdo.
Other minor corrections for HDO cycle.
Transition from fractionation coefficients from Merlivat et al. 1967 to Lamb et al. 2017

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     &                      qperemin
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.qperemin) 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?
107c               alpha_c(ig) = exp(16288./zt(ig,1)**2.-9.34e-2)
108                alpha_c = exp(13525./zt(ig,1)**2.-5.59e-2) !Lamb
109               else
110                alpha_c(ig) = 1.
111               endif
112
113               if (old_h2o_vap(ig).gt.qperemin) 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.