Changeset 2312 for trunk/LMDZ.MARS/libf/phymars/vdifc_mod.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/vdifc_mod.F (modified) (14 diffs)

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