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

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

File:

• trunk/LMDZ.MARS/libf/phymars/simpleclouds.F (modified) (6 diffs)

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