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

Last change on this file since 2461 was 2378, checked in by lrossi, 5 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
RevLine 
[2312]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,
[2316]9     &                      old_h2o_vap,pdqsdif,dwatercap_dif,
[2312]10     &                      hdoflux)
11
12      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice,
[2324]13     &                      igcm_hdo_vap, igcm_hdo_ice,
14     &                      qperemin
[2312]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)
[2316]39      REAL, INTENT(IN) :: dwatercap_dif(ngrid)  ! trend related to permanent ice
[2312]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
[2316]52      REAL h2oflux(ngrid)       ! value of vapour flux of H2O
53                                      ! same sign as pdqsdif
[2312]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
[2316]65             
66            h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice) +
67     &          dwatercap_dif(ig)
[2312]68
69            !! IF Sublimation
70            if (h2oflux(ig).le.0.) then
71
[2324]72               if (pqsurf(ig,igcm_h2o_ice).gt.qperemin) then
[2312]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
[2316]93C               dwatercap_dif is how much we sublimate in excess of
[2312]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
[2316]98c               dwatercap_dif* 5 * 2 * 155.76e-6 (=1 SMOW)
[2312]99                    hdoflux(ig)= hdoflux(ig)
[2316]100     &                   +(dwatercap_dif(ig)*(2.*155.76e-6)*5.)
101                endif
102             endif ! watercap
[2312]103
104            else ! condensation
105
106               if (hdofrac) then !do we use fractionation?
[2378]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
[2312]109               else
110                alpha_c(ig) = 1.
111               endif
112
[2324]113               if (old_h2o_vap(ig).gt.qperemin) then
[2312]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.