Changeset 2312 for trunk/LMDZ.MARS/libf/phymars

Timestamp:

May 6, 2020, 4:46:33 PM (5 years ago)

Author:

lrossi

Message:

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

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• callkeys.h (modified) (2 diffs)
• callsedim_mod.F (modified) (2 diffs)
• conf_phys.F (modified) (2 diffs)
• dyn1d/testphys1d.F (modified) (4 diffs)
• hdo_surfex_mod.F (added)
• initracer.F (modified) (7 diffs)
• phyetat0_mod.F90 (modified) (2 diffs)
• phyredem.F90 (modified) (3 diffs)
• physiq_mod.F (modified) (25 diffs)
• simpleclouds.F (modified) (6 diffs)
• tracer_mod.F90 (modified) (1 diff)
• vdifc_mod.F (modified) (14 diffs)
• watercloud_mod.F (modified) (9 diffs)

trunk/LMDZ.MARS/libf/phymars/callkeys.h

@@ -16 +16 @@
      &   ,co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying        &
      &   ,satindexco2,rdstorm,slpwind,calllott_nonoro                   &
-     &   ,latentheat_surfwater,gwd_convective_source,startphy_file
-
+     &   ,latentheat_surfwater,gwd_convective_source,startphy_file      &
+     &   ,hdo,hdofrac
+     
       COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
      &   ,dustbin,nltemodel,nircorr,solvarmod,solvaryear,dustinjection
@@ -78 +79 @@
       logical CLFvarying
       logical water
+      logical hdo
+      logical hdofrac
       logical microphys
       logical photochem
       integer nltemodel
       integer nircorr
-
 
       character(len=100) dustiropacity

trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F

@@ -17 +17 @@
      &                      rho_dust, rho_q, radius, varian,
      &                      igcm_ccn_mass, igcm_ccn_number,
-     &                      igcm_h2o_ice, nuice_sed, nuice_ref,
+     &                      igcm_h2o_ice, igcm_hdo_ice,
+     &                      nuice_sed, nuice_ref,
      &                      igcm_ccnco2_mass,igcm_ccnco2_number,
      &                      igcm_co2_ice, igcm_stormdust_mass, 
@@ -499 +500 @@
 c -----------------------------------------------------------------
            else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
-     &       .or. (iq .eq. igcm_h2o_ice)) then
+     &       .or. (iq .eq. igcm_h2o_ice)
+     &       .or. (iq .eq. igcm_hdo_ice)) then
             if (microphys) then
               ! water ice sedimentation

trunk/LMDZ.MARS/libf/phymars/conf_phys.F

@@ -541 +541 @@
          call getin_p("water",water)
          write(*,*) " water = ",water
+! hdo
+         write(*,*) "Compute hdo cycle ?"
+         hdo=.false. ! default value
+         call getin_p("hdo",hdo)
+         write(*,*) " hdo = ",hdo
+
+         write(*,*) "Use fractionation for hdo?"
+         hdofrac=.true. ! default value
+         call getin_p("hdofrac",hdofrac)
+         write(*,*) " hdofrac = ",hdofrac
+
 
 ! sub-grid cloud fraction: fixed
@@ -618 +629 @@
      &          "water requires tracer",1)
          endif
+
+         if (hdo.and..not.water) then
+           print*,'if hdo is used, water should be used too'
+           stop
+         endif
+
          
 ! water ice clouds effective variance distribution for sedimentaion       

trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F

@@ -117 +117 @@
       integer :: nq=1 ! number of tracers
       real :: latitude(1), longitude(1), cell_area(1)
+      integer iqh2ovap
+      integer iqh2oice
 
       character*2 str2
@@ -294 +296 @@
           ! WATER VAPOUR
           if (txt.eq."h2o_vap") then
+             iqh2ovap=iq !remember index for water vap
             !look for a "profile_h2o_vap" input file
             open(91,file='profile_h2o_vap',status='old',
@@ -309 +312 @@
           ! WATER ICE
           if (txt.eq."h2o_ice") then
+                iqh2oice = iq !remember index for water ice
             !look for a "profile_h2o_vap" input file
             open(91,file='profile_h2o_ice',status='old',
@@ -322 +326 @@
             close(91)
           endif ! of if (txt.eq."h2o_ice")
+
+          ! HDO VAPOUR
+          if (txt.eq."hdo_vap") then
+            !look for a "profile_hdo_vap" input file
+            open(91,file='profile_hdo_vap',status='old',
+     &       form='formatted',iostat=ierr)
+            if (ierr.eq.0) then
+              read(91,*) qsurf(iq)
+              do ilayer=1,nlayer
+                read(91,*) q(ilayer,iq)
+              enddo
+            else
+              write(*,*) "No profile_hdo_vap file!"
+              do ilayer=1,nlayer
+                q(ilayer,iq) =  q(ilayer,iqh2ovap)*2*155.76e-6*5 
+              enddo
+            endif
+            close(91)
+          endif ! of if (txt.eq."hdo_ice")
+          ! HDO ICE
+          if (txt.eq."hdo_ice") then
+            !look for a "profile_hdo_vap" input file
+            open(91,file='profile_hdo_ice',status='old',
+     &       form='formatted',iostat=ierr)
+            if (ierr.eq.0) then
+              read(91,*) qsurf(iq)
+              do ilayer=1,nlayer
+                read(91,*) q(ilayer,iq)
+              enddo
+            else
+              write(*,*) "No profile_hdo_ice file!"
+                qsurf(iq) = qsurf(iqh2oice) * 2*155.76e-6*5
+               do ilayer=1,nlayer
+                q(ilayer,iq) = q(ilayer,iqh2oice) * 2*155.76e-6*5 
+              enddo
+            endif
+            close(91)
+          endif ! of if (txt.eq."hdo_ice")
+
+
           ! DUST
           !if (txt(1:4).eq."dust") then

trunk/LMDZ.MARS/libf/phymars/initracer.F

@@ -70 +70 @@
       igcm_h2o_vap=0
       igcm_h2o_ice=0
+      igcm_hdo_vap=0
+      igcm_hdo_ice=0
       igcm_stormdust_mass=0
       igcm_stormdust_number=0
@@ -345 +347 @@
           count=count+1
         endif
+        if (noms(iq).eq."hdo_vap") then
+          igcm_hdo_vap=iq
+          mmol(igcm_hdo_vap)=18. !19
+          count=count+1
+        endif
         if (noms(iq).eq."co2_ice") then
           igcm_co2_ice=iq
@@ -353 +360 @@
           igcm_h2o_ice=iq
           mmol(igcm_h2o_ice)=18.
+          count=count+1
+        endif
+        if (noms(iq).eq."hdo_ice") then
+          igcm_hdo_ice=iq
+          mmol(igcm_hdo_ice)=18. !19
           count=count+1
         endif
@@ -403 +415 @@
       endif
 
+      ! Additional test required for HDO
+      ! We need to compute some things for H2O before HDO
+      if (hdo) then
+        if (igcm_h2o_vap.gt.igcm_hdo_vap) then
+           write(*,*) "Tracer H2O must be initialized before HDO"
+           STOP
+        endif
+      endif
+
 c------------------------------------------------------------
 c     Initialize tracers .... (in tracer_mod)
@@ -429 +450 @@
 c       iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
 
-        if( (nq.lt.2).or.(water.and.(nq.lt.4)) ) then
+        if( (nq.lt.2).or.(water.and.(nq.lt.4))
+     *       .or.(hdo.and.(nq.lt.6) )) then
           write(*,*)'initracer: nq is too low : nq=', nq
           write(*,*)'water= ',water,' doubleq= ',doubleq   
@@ -594 +616 @@
 
       end if  ! (water)
+
+c     Initialization for hdo vapor
+c     ------------------------------
+      if (hdo) then
+         radius(igcm_hdo_vap)=0.
+         alpha_lift(igcm_hdo_vap) =0.
+         alpha_devil(igcm_hdo_vap)=0.
+         if(water.and.(nq.ge.2)) then
+           radius(igcm_hdo_ice)=3.e-6
+           rho_q(igcm_hdo_ice)=rho_ice
+           alpha_lift(igcm_hdo_ice) =0.
+           alpha_devil(igcm_hdo_ice)=0.
+         elseif(hdo.and.(nq.lt.6)) then
+            write(*,*) 'nq is too low : nq=', nq
+            write(*,*) 'hdo= ',hdo
+         endif
+
+      end if  ! (hdo)
+
       
 ! Initialisation for CO2 clouds
@@ -683 +724 @@
          endif
        endif
+
+       if (hdo) then
+       ! verify that we indeed have hdo_vap and hdo_ice tracers
+         if (igcm_hdo_vap.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use hdo option without ",
+     &                "an hdo_vap tracer !"
+           stop
+         endif
+         if (igcm_hdo_ice.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use hdo option without ",
+     &                "an hdo_ice tracer !"
+           stop
+         endif
+       endif
+
 
        if (co2clouds) then

trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90

@@ -293 +293 @@
             minval(base), maxval(base)
 
- 
 ! Time axis
 ! obtain timestart from run.def
@@ -567 +566 @@
       write(*,*) 'iq = ', iq
     endif
+
+    if (hdo) then
+    if (txt.eq."hdo_vap") then
+      txt="hdo_ice"
+      write(*,*) 'phyetat0: loading surface tracer', &
+                           ' hdo_ice instead of hdo_vap'
+    endif
+    endif !hdo
+
     if (startphy_file) then
        call get_field(txt,qsurf(:,iq),found,indextime)

trunk/LMDZ.MARS/libf/phymars/phyredem.F90

@@ -189 +189 @@
   integer :: i_h2o_vap=0
   integer :: i_h2o_ice=0
+  integer :: i_hdo_vap=0
+  integer :: i_hdo_ice=0
 
   
@@ -245 +247 @@
       i_h2o_ice=iq
     endif
+
+  ! look for HDO vapour & HDO ice tracers (if any)
+    if (noms(iq).eq."hdo_vap") then
+      i_hdo_vap=iq
+    endif
+    if (noms(iq).eq."hdo_ice") then
+      i_hdo_ice=iq
+    endif
   enddo
+
   
   if (nq.gt.0) then
@@ -264 +275 @@
         endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
       endif ! of if (txt.eq."h2o_vap")
+
+      if (txt.eq."hdo_vap") then
+        write(*,*)"physdem1: skipping HDO vapour tracer"
+        if (i_hdo_ice.eq.i_hdo_vap) then
+          ! then there is no "water ice" tracer; but still
+          ! there is some water ice on the surface
+          write(*,*)"          writing HDO ice instead"
+          txt="hdo_ice"
+        else
+          ! there is a "water ice" tracer which has been / will be
+          ! delt with in due time
+          cycle
+        endif ! of if (igcm_hdo_ice.eq.igcm_hdo_vap)
+      endif ! of if (txt.eq."hdo_vap")
+
       call put_field(trim(txt),"tracer on surface",qsurf(:,iq),time)
     enddo

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -27 +27 @@
       use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice,
      &                      igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
+     &                      igcm_hdo_vap, igcm_hdo_ice,
      &                      igcm_ccn_mass, igcm_ccn_number,
      &                      igcm_ccnco2_mass, igcm_ccnco2_number,
@@ -87 +88 @@
       use wxios, only: wxios_context_init, xios_context_finalize
 #endif
+
 
       IMPLICIT NONE
@@ -395 +397 @@
                                !   Qabs instead of Qext
                                !   (direct comparison with TES)
+      REAL mtotD(ngrid)          ! Total mass of HDO vapor (kg/m2)
+      REAL icetotD(ngrid)        ! Total mass of HDO ice (kg/m2)
+      REAL DoH_vap(ngrid,nlayer) !D/H ratio
+      REAL DoH_ice(ngrid,nlayer) !D/H ratio
+      REAL DoH_surf(ngrid) !D/H ratio surface
                                
       REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
@@ -418 +425 @@
 c Test 1d/3d scavenging
       real h2otot(ngrid)
+      real hdotot(ngrid)
       REAL satu(ngrid,nlayer)  ! satu ratio for output
       REAL zqsat(ngrid,nlayer) ! saturation
@@ -701 +709 @@
       call update_xios_timestep
 #endif      
+
+c     Initialize various variables
 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       pdv(:,:)=0
@@ -714 +724 @@
       dsotop(:,:)=0.
       dwatercap(:)=0
-
+      
 #ifdef DUSTSTORM
       pq_tmp(:,:,:)=0
@@ -773 +783 @@
          ENDDO
       ENDDO
+
 #ifndef MESOSCALE
 c-----------------------------------------------------------------------
@@ -886 +897 @@
            !! callradite for background dust in the case of slope wind entrainment
            nohmons=.true.
+
 c          Radiative transfer
 c          ------------------
@@ -1200 +1212 @@
            endif ! end of if (dustinjection.gt.0)
 
+
 c-----------------------------------------------------------------------
 c    4. Gravity wave and subgrid scale topography drag :
@@ -1263 +1276 @@
              pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
          ENDIF
+
 c        Calling vdif (Martian version WITH CO2 condensation)
          CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
@@ -1289 +1303 @@
             ENDDO
           ENDDO
+
           if (tracer) then
            DO iq=1, nq
@@ -1369 +1384 @@
           q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
         ENDDO      
-    
+
         if (tracer) then
         DO iq=1,nq
@@ -1388 +1403 @@
         lmax_th_out(:)=0.
       end if
+
 c-----------------------------------------------------------------------
 c   7. Dry convective adjustment:
@@ -1483 +1499 @@
      &          zdtcloud(1:ngrid,1:nlayer)
            endif
-           
+
 ! increment water vapour and ice atmospheric tracers tendencies
            pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = 
@@ -1491 +1507 @@
      &       pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + 
      &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
+
+           if (hdo) then
+! increment HDO vapour and ice atmospheric tracers tendencies
+           pdq(1:ngrid,1:nlayer,igcm_hdo_vap) = 
+     &       pdq(1:ngrid,1:nlayer,igcm_hdo_vap) + 
+     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap)
+           pdq(1:ngrid,1:nlayer,igcm_hdo_ice) = 
+     &       pdq(1:ngrid,1:nlayer,igcm_hdo_ice) + 
+     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice)
+           endif !hdo
 
 ! increment dust and ccn masses and numbers
@@ -1896 +1922 @@
             ENDDO
            ENDDO
-                   
-           DO iq=1, nq
+
+          DO iq=1, nq
             DO ig=1,ngrid
               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqsc(ig,iq)
             ENDDO  ! (ig)
            ENDDO    ! (iq)
-           
+
          ENDIF ! of IF (tracer)
 
@@ -2248 +2274 @@
      &                dsotop(ig,:)*tauscaling(ig)
               enddo
+         
            if(doubleq) then
               do ig=1,ngrid 
@@ -2325 +2352 @@
              rave(:)=0
              tauTES(:)=0
+
+             IF (hdo) then
+                 mtotD(:)=0
+                 icetotD(:)=0
+             ENDIF !hdo
+
              do ig=1,ngrid 
                do l=1,nlayer
@@ -2333 +2366 @@
      &                        zq(ig,l,igcm_h2o_ice) * 
      &                        (zplev(ig,l) - zplev(ig,l+1)) / g
+                 IF (hdo) then
+                 mtotD(ig) = mtotD(ig) + 
+     &                      zq(ig,l,igcm_hdo_vap) * 
+     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
+                 icetotD(ig) = icetotD(ig) + 
+     &                        zq(ig,l,igcm_hdo_ice) * 
+     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
+                 ENDIF !hdo
+
 c                Computing abs optical depth at 825 cm-1 in each
 c                  layer to simulate NEW TES retrieval
@@ -2887 +2929 @@
      &             3,zq(:,:,igcm_h2o_vap))
 
+            if (hdo) then
+            vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice)
+     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice)
+            CALL WRITEDIAGFI(ngrid,'vmr_hdoice','hdo ice vmr',
+     &                       'mol/mol',3,vmr)
+            vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap)
+     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap)
+            CALL WRITEDIAGFI(ngrid,'vmr_hdovap','hdo vap vmr',
+     &                       'mol/mol',3,vmr)
+            call WRITEDIAGFI(ngrid,'hdo_ice','hdo_ice','kg/kg',
+     &             3,zq(:,:,igcm_hdo_ice))
+            call WRITEDIAGFI(ngrid,'hdo_vap','hdo_vap','kg/kg',
+     &             3,zq(:,:,igcm_hdo_vap))
+
+            CALL WRITEDIAGFI(ngrid,'mtotD',
+     &                       'total mass of HDO vapor',
+     &                       'kg/m2',2,mtotD)
+            CALL WRITEDIAGFI(ngrid,'icetotD',
+     &                       'total mass of HDO ice',
+     &                       'kg/m2',2,icetotD)
+
+C           Calculation of the D/H ratio
+            do l=1,nlayer
+                do ig=1,ngrid
+                if (zq(ig,l,igcm_h2o_vap).gt.1e-16) then 
+                    DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
+     &              zq(ig,l,igcm_h2o_vap) )*1./(2.*155.76e-6)
+                else
+                    DoH_vap(ig,l) = 0.
+                endif
+                enddo
+            enddo
+
+            do l=1,nlayer
+                do ig=1,ngrid
+                if (zq(ig,l,igcm_h2o_ice).gt.1e-16) then
+                    DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/
+     &                  zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6)
+                else
+                    DoH_ice(ig,l) = 0.
+                endif
+                enddo
+            enddo
+
+            CALL WRITEDIAGFI(ngrid,'DoH_vap',
+     &                       'D/H ratio in vapor',
+     &                       ' ',3,DoH_vap)  
+            CALL WRITEDIAGFI(ngrid,'DoH_ice',
+     &                       'D/H ratio in ice',
+     &                       '',3,DoH_ice)
+
+            endif !hdo
 
 !A. Pottier
@@ -2907 +3001 @@
      &                       'surface h2o_ice',
      &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
+            if (hdo) then
+            call WRITEDIAGFI(ngrid,'hdo_ice_s',
+     &                       'surface hdo_ice',
+     &                       'kg.m-2',2,qsurf(1,igcm_hdo_ice))
+
+                do ig=1,ngrid
+                if (qsurf(ig,igcm_h2o_ice).gt.1e-30) then
+                    DoH_surf(ig) = 0.5*( qsurf(ig,igcm_hdo_ice)/
+     &                  qsurf(ig,igcm_h2o_ice) )/155.76e-6
+                else
+                    DoH_surf(ig) = 0.
+                endif
+                enddo
+
+            call WRITEDIAGFI(ngrid,'DoH_surf',
+     &                       'surface D/H',
+     &                       '',2,DoH_surf)
+            endif ! hdo
+
             CALL WRITEDIAGFI(ngrid,'albedo',
      &                         'albedo',
@@ -3430 +3543 @@
            icetot = 0
            h2otot = qsurf(1,igcm_h2o_ice)
+           if (hdo) THEN
+           mtotD = 0
+           icetotD = 0
+           hdotot = qsurf(1,igcm_hdo_ice)
+           ENDIF !hdo
 
            do l=1,nlayer
@@ -3436 +3554 @@
              icetot = icetot +  zq(1,l,igcm_h2o_ice) 
      &                 * (zplev(1,l) - zplev(1,l+1)) / g
+               if (hdo) THEN
+                 mtotD = mtotD +  zq(1,l,igcm_hdo_vap) 
+     &                 * (zplev(1,l) - zplev(1,l+1)) / g
+                 icetotD = icetotD +  zq(1,l,igcm_hdo_ice) 
+     &                 * (zplev(1,l) - zplev(1,l+1)) / g
+               ENDIF !hdo
            end do
            h2otot = h2otot+mtot+icetot
+               IF (hdo) then
+                   hdotot = hdotot+mtotD+icetotD
+               ENDIF ! hdo
+                
 
              CALL WRITEDIAGFI(ngrid,'h2otot',
@@ -3448 +3576 @@
      &                         'icetot',
      &                         'kg/m2',0,icetot)
+
+             IF (hdo) THEN
+             CALL WRITEDIAGFI(ngrid,'mtotD',
+     &                         'mtotD',
+     &                         'kg/m2',0,mtotD)
+             CALL WRITEDIAGFI(ngrid,'icetotD',
+     &                         'icetotD',
+     &                         'kg/m2',0,icetotD)
+             CALL WRITEDIAGFI(ngrid,'hdotot',
+     &                         'hdotot',
+     &                         'kg/m2',0,hdotot)
+
+C           Calculation of the D/H ratio
+            do l=1,nlayer
+                if (zq(1,l,igcm_h2o_vap).gt.1e-16) then
+                    DoH_vap(1,l) = 0.5*( zq(1,l,igcm_hdo_vap)/
+     &                zq(1,l,igcm_h2o_vap) )/155.76e-6
+                else
+                    DoH_vap(1,l) = 0.
+                endif
+            enddo
+
+            do l=1,nlayer
+                if (zq(1,l,igcm_h2o_ice).gt.1e-16) then
+                    DoH_ice(1,l) = 0.5*( zq(1,l,igcm_hdo_ice)/
+     &                  zq(1,l,igcm_h2o_ice) )/155.76e-6
+                else
+                    DoH_ice(1,l) = 0.
+                endif
+            enddo
+
+            CALL WRITEDIAGFI(ngrid,'DoH_vap',
+     &                       'D/H ratio in vapor',
+     &                       ' ',1,DoH_vap)  
+            CALL WRITEDIAGFI(ngrid,'DoH_ice',
+     &                       'D/H ratio in ice',
+     &                       '',1,DoH_ice)
+
+             ENDIF !Hdo
+
 
            if (scavenging) then
@@ -3527 +3695 @@
      &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
      &                          +zdqcloud(1,:,igcm_h2o_vap))
+        IF (hdo) THEN
+        call WRITEDIAGFI(ngrid,'zdqcloud_iceD','cloud ice hdo',
+     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_ice))
+        call WRITEDIAGFI(ngrid,'zdqcloud_vapD','cloud vap hdo',
+     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_vap))
+
+        ENDIF ! hdo
      
         call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,

trunk/LMDZ.MARS/libf/phymars/simpleclouds.F

@@ -5 +5 @@
       USE updaterad
       USE watersat_mod, ONLY: watersat
-      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice
+      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice,
+     &                      igcm_hdo_vap, igcm_hdo_ice
       USE comcstfi_h
       use dimradmars_mod, only: naerkind
+
       implicit none
 c------------------------------------------------------------------
@@ -77 +79 @@
       REAL ccntyp(ngrid,nlay)
                                       ! Typical dust number density (#/kg)
+      REAL alpha_c(ngrid,nlay) !!MARGAUX: alpha_c as a spatial variable
+
 c     CCN reduction factor
 c      REAL, PARAMETER :: ccn_factor = 4.5  !! comme TESTS_JB // 1. avant
-      
-
+      REAL DoH_vap(ngrid,nlay)
+      REAL DoH_ice(ngrid,nlay)
 c-----------------------------------------------------------------------
 c    1. initialisation
@@ -100 +104 @@
           zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004 
           zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
+
+          if (hdo) then
+          zq(ig,l,igcm_hdo_vap)=
+     &      pq(ig,l,igcm_hdo_vap)+pdq(ig,l,igcm_hdo_vap)*ptimestep
+          zq(ig,l,igcm_hdo_vap)=max(zq(ig,l,igcm_hdo_vap),1e-30) ! FF 12/2004 
+          zq0(ig,l,igcm_hdo_vap)=zq(ig,l,igcm_hdo_vap)
+
+          zq(ig,l,igcm_hdo_ice)=
+     &      pq(ig,l,igcm_hdo_ice)+pdq(ig,l,igcm_hdo_ice)*ptimestep
+          zq(ig,l,igcm_hdo_ice)=max(zq(ig,l,igcm_hdo_ice),1e-30) ! FF 12/2004 
+          zq0(ig,l,igcm_hdo_ice)=zq(ig,l,igcm_hdo_ice)
+
+          endif !hdo
         enddo
       enddo
 
-
       pdqcloud(1:ngrid,1:nlay,1:nq)=0
       pdtcloud(1:ngrid,1:nlay)=0
 
+      alpha_c(1:ngrid,1:nlay)=0.
+
 c     ----------------------------------------------
-c
 c
 c     Rapport de melange a saturation dans la couche l : -------
@@ -122 +139 @@
 
           if (zq(ig,l,igcm_h2o_vap).ge.zqsat(ig,l))then  !  Condensation
-            dzq=zq(ig,l,igcm_h2o_vap)-zqsat(ig,l)                
+            dzq=zq(ig,l,igcm_h2o_vap)-zqsat(ig,l)
+
           elseif(zq(ig,l,igcm_h2o_vap).lt.zqsat(ig,l))then  ! Sublimation
             dzq=-min(zqsat(ig,l)-zq(ig,l,igcm_h2o_vap),
      &               zq(ig,l,igcm_h2o_ice))
           endif
-
+            
 c         Water Mass change
 c         ~~~~~~~~~~~~~~~~~
           zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)+dzq
           zq(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)-dzq
-          
-
+        
         enddo ! of do ig=1,ngrid
       enddo ! of do l=1,nlay
+
 
 c     Tendance finale
@@ -145 +163 @@
           pdqcloud(ig,l,igcm_h2o_ice) =
      &      (zq(ig,l,igcm_h2o_ice) - zq0(ig,l,igcm_h2o_ice))/ptimestep
+
           lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
           pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
-        end do
-      end do
+
+          if (hdo) then
+            if (pdqcloud(ig,l,igcm_h2o_ice).gt.0.0) then !condens
+
+                if (hdofrac) then ! do we use fractionation?
+                alpha_c(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
+c               alpha_c(ig,l) = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
+                else
+                alpha_c(ig,l) = 1.d0
+                endif
+               
+                if (zq0(ig,l,igcm_h2o_vap).gt.1.e-16) then
+                  pdqcloud(ig,l,igcm_hdo_ice)= 
+     &              pdqcloud(ig,l,igcm_h2o_ice)*alpha_c(ig,l)*
+     &         ( zq0(ig,l,igcm_hdo_vap)
+     &                 /zq0(ig,l,igcm_h2o_vap) )
+                else
+                   pdqcloud(ig,l,igcm_hdo_ice)= 0.0
+                endif
+
+                pdqcloud(ig,l,igcm_hdo_ice) =
+     &            min(pdqcloud(ig,l,igcm_hdo_ice),
+     &              zq0(ig,l,igcm_hdo_vap)/ptimestep)
+
+                pdqcloud(ig,l,igcm_hdo_vap)= 
+     &               -pdqcloud(ig,l,igcm_hdo_ice)       
+
+          else  ! sublimation
+
+             if (zq0(ig,l,igcm_h2o_ice).gt.1.e-16) then 
+                pdqcloud(ig,l,igcm_hdo_ice)= 
+     &               pdqcloud(ig,l,igcm_h2o_ice)*
+     &      ( zq0(ig,l,igcm_hdo_ice)
+     &              /zq0(ig,l,igcm_h2o_ice) )
+             else
+                pdqcloud(ig,l,igcm_hdo_ice)= 0.
+             endif
+
+              pdqcloud(ig,l,igcm_hdo_ice) =
+     &          max(pdqcloud(ig,l,igcm_hdo_ice),
+     &            -zq0(ig,l,igcm_hdo_ice)/ptimestep)
+
+              pdqcloud(ig,l,igcm_hdo_vap)= 
+     &             -pdqcloud(ig,l,igcm_hdo_ice)        
+
+            endif ! condensation/sublimation
+
+          endif ! hdo
+
+        enddo ! of do ig=1,ngrid
+      enddo ! of do l=1,nlay
 
 c     ice crystal radius
@@ -158 +226 @@
       end do
 
+c     if (hdo) then
+c           CALL WRITEDIAGFI(ngrid,'alpha_c',
+c    &                       'alpha_c',
+c    &                       ' ',3,alpha_c) 
+c     endif !hdo
 c------------------------------------------------------------------
       return

trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90

@@ -53 +53 @@
       integer,save :: igcm_h2o_vap ! water vapour
       integer,save :: igcm_h2o_ice ! water ice
+      integer,save :: igcm_hdo_vap ! hdo vapour
+      integer,save :: igcm_hdo_ice ! hdo ice
       integer,save :: igcm_co2_ice ! co2 ice
 

trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F

@@ -18 +18 @@
      &                      igcm_dust_submicron, igcm_h2o_vap,
      &                      igcm_h2o_ice, alpha_lift,
+     &                      igcm_hdo_vap, igcm_hdo_ice,
      &                      igcm_stormdust_mass, igcm_stormdust_number
       use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness
@@ -24 +25 @@
       use turb_mod, only: turb_resolved, ustar, tstar
       use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol
+      use hdo_surfex_mod, only: hdo_surfex
+
       IMPLICIT NONE
 
@@ -140 +143 @@
       REAL qsat(ngrid)
 
+      REAL hdoflux(ngrid)       ! value of vapour flux of HDO
+      REAL h2oflux(ngrid)       ! value of vapour flux of H2O
+      REAL old_h2o_vap(ngrid)   ! traceur d'eau avant traitement
+
       REAL kmixmin
 
@@ -167 +174 @@
 
       REAL,INTENT(OUT) :: sensibFlux(ngrid)
+
+!!MARGAUX
+      REAL DoH_vap(ngrid,nlay)
 
 c    ** un petit test de coherence
@@ -443 +453 @@
 c       donc les entrees sont /zcdv/ pour la condition a la limite sol
 c       et /zkv/ = Ku
- 
+
       zb(1:ngrid,2:nlay)=zkv(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
       zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
@@ -786 +796 @@
           zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
 
-          if ((water).and.(iq.eq.igcm_h2o_vap)) then 
+          if ((water).and.(iq.eq.igcm_h2o_vap)) then
 c            This line is required to account for turbulent transport 
 c            from surface (e.g. ice) to mid-layer of atmosphere:
@@ -811 +821 @@
           ENDDO
 
-          if (water.and.(iq.eq.igcm_h2o_ice)) then
+          if ( (water.and.(iq.eq.igcm_h2o_ice)) 
+     $      .or. (hdo.and.(iq.eq.igcm_hdo_ice)) ) then
             ! special case for water ice tracer: do not include
             ! h2o ice tracer from surface (which is set when handling
@@ -821 +832 @@
      $         zb(ig,2)*zc(ig,2)) *z1(ig)
             ENDDO
+         else if (hdo.and.(iq.eq.igcm_hdo_vap)) then
+            CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
+     &                      zt,zq,pqsurf,
+     &                      old_h2o_vap,pdqsdif,h2oflux,
+     &                      hdoflux)
+            DO ig=1,ngrid
+                z1(ig)=1./(za(ig,1)+zb(ig,1)+
+     $           zb(ig,2)*(1.-zd(ig,2)))
+                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
+     $         zb(ig,2)*zc(ig,2) +
+     $        (-hdoflux(ig)) *ptimestep) *z1(ig)  !tracer flux from surface
+            ENDDO
+
           else ! general case
             DO ig=1,ngrid
@@ -829 +853 @@
      $        (-pdqsdif(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
             ENDDO
+
           endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
   
@@ -834 +859 @@
 c           Calculation for turbulent exchange with the surface (for ice)
             DO ig=1,ngrid
+              old_h2o_vap(ig)=zq(ig,1,igcm_h2o_vap)
               zd(ig,1)=zb(ig,1)*z1(ig)
               zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
@@ -843 +869 @@
 
             DO ig=1,ngrid
+            IF (hdo) then !if hdo, need to treat polar caps differently
+              h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice) ! h2oflux is
+                                                     ! uncorrected waterflux
+            ENDIF !hdo
+
               if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
      &             .gt.(pqsurf(ig,igcm_h2o_ice))) then
@@ -857 +888 @@
      &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
                    dwatercap_dif(ig) = 0.0
+                   if (hdo) then
+                     h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice)
+                   endif
 
                 endif   
@@ -910 +944 @@
         enddo ! of do iq=1,nq
       end if ! of if(tracer)
-      
+
+C       Diagnostic output for HDO
+        if (hdo) then
+            CALL WRITEDIAGFI(ngrid,'hdoflux',
+     &                       'hdoflux',
+     &                       ' ',2,hdoflux)  
+            CALL WRITEDIAGFI(ngrid,'h2oflux',
+     &                       'h2oflux',
+     &                       ' ',2,h2oflux)
+        endif
 
 c-----------------------------------------------------------------------
@@ -958 +1001 @@
       END SUBROUTINE vdifc
 
+c====================================
+
+
       END MODULE vdifc_mod

trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F

@@ -19 +19 @@
       USE watersat_mod, ONLY: watersat
       use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
+     &                      igcm_hdo_vap, igcm_hdo_ice,
      &                      igcm_dust_mass, igcm_dust_number,
      &                      igcm_ccn_mass, igcm_ccn_number,
@@ -124 +125 @@
       REAL :: mincloud ! min cloud frac
       LOGICAL, save :: flagcloud=.true. 
+
 c    ** un petit test de coherence
 c       --------------------------
@@ -312 +314 @@
      &          sum_subpdq(ig,l,igcm_h2o_vap)
      &        + pdq(ig,l,igcm_h2o_vap)
+            IF (hdo) THEN
+            sum_subpdq(ig,l,igcm_hdo_ice) = 
+     &          sum_subpdq(ig,l,igcm_hdo_ice)
+     &        + pdq(ig,l,igcm_hdo_ice)
+            sum_subpdq(ig,l,igcm_hdo_vap) = 
+     &          sum_subpdq(ig,l,igcm_hdo_vap)
+     &        + pdq(ig,l,igcm_hdo_vap)
+            ENDIF !hdo
           ENDDO
         ENDDO      
@@ -361 +371 @@
      &          sum_subpdq(ig,l,igcm_h2o_vap)
      &        + subpdqcloud(ig,l,igcm_h2o_vap)
+
+            IF (hdo) THEN
+            sum_subpdq(ig,l,igcm_hdo_ice) =
+     &          sum_subpdq(ig,l,igcm_hdo_ice)
+     &        + subpdqcloud(ig,l,igcm_hdo_ice)
+            sum_subpdq(ig,l,igcm_hdo_vap) =
+     &          sum_subpdq(ig,l,igcm_hdo_vap)
+     &        + subpdqcloud(ig,l,igcm_hdo_vap)
+            ENDIF ! hdo
+
           ENDDO
         ENDDO
@@ -397 +417 @@
      &        sum_subpdq(ig,l,igcm_h2o_vap)/real(imicro) 
      &       - pdq(ig,l,igcm_h2o_vap)
+            IF (hdo) THEN
+            pdqcloud(ig,l,igcm_hdo_ice) = 
+     &        sum_subpdq(ig,l,igcm_hdo_ice)/real(imicro) 
+     &       - pdq(ig,l,igcm_hdo_ice)
+            pdqcloud(ig,l,igcm_hdo_vap) = 
+     &        sum_subpdq(ig,l,igcm_hdo_vap)/real(imicro) 
+     &       - pdq(ig,l,igcm_hdo_vap)
+            ENDIF !hdo
          ENDDO
        ENDDO
@@ -478 +506 @@
         DO l=1,nlay
          DO ig=1,ngrid
+
           IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
      &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice)) 
@@ -484 +513 @@
      &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
            pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
+             if (hdo) then
+           pdqcloud(ig,l,igcm_hdo_ice) = 
+     &     - pq(ig,l,igcm_hdo_ice)/ptimestep - pdq(ig,l,igcm_hdo_ice)
+           pdqcloud(ig,l,igcm_hdo_vap) = -pdqcloud(ig,l,igcm_hdo_ice)
+             endif
           ENDIF 
           IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
@@ -491 +525 @@
      &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
            pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
+             if (hdo) then
+           pdqcloud(ig,l,igcm_hdo_vap) = 
+     &     - pq(ig,l,igcm_hdo_vap)/ptimestep - pdq(ig,l,igcm_hdo_vap)
+           pdqcloud(ig,l,igcm_hdo_ice) = -pdqcloud(ig,l,igcm_hdo_vap)
+             endif
           ENDIF
-         ENDDO
-        ENDDO
-
+
+         ENDDO
+        ENDDO
 
 c------Update the ice and dust particle size "rice" for output or photochemistry
@@ -627 +666 @@
      &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
      &      igcm_h2o_vap))
+      if (hdo) then
+      call WRITEDIAGFI(ngrid,"pdqiceD","pdqiceD apres microphysique"
+     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,igcm_hdo_ice))
+      call WRITEDIAGFI(ngrid,"pdqvapD","pdqvapD apres microphysique"
+     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
+     &      igcm_hdo_vap))
+      endif
       call WRITEDIAGFI(ngrid,"pdqccn2","pdqcloudccn apres microphysique"
      &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,