Changeset 2316 for trunk/LMDZ.MARS/libf

Timestamp:

May 12, 2020, 11:51:22 AM (5 years ago)

Author:

lrossi

Message:

Mars GCM:
Fixing some errors in vdifc_mod related to variable watercap. This variable was also integrated to the hdo cycle.
Also added watercap output for the 1D model.
LR

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• hdo_surfex_mod.F (modified) (5 diffs)
• physiq_mod.F (modified) (1 diff)
• vdifc_mod.F (modified) (4 diffs)

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

@@ -7 +7 @@
       subroutine hdo_surfex(ngrid,nlay,nq,ptimestep,
      &                      zt,zq,pqsurf,
-     &                      old_h2o_vap,pdqsdif,h2oflux,
+     &                      old_h2o_vap,pdqsdif,dwatercap_dif,
      &                      hdoflux)
 
@@ -37 +37 @@
       REAL, INTENT(IN) :: old_h2o_vap(ngrid)     ! traceur d'eau avant
                                            !traitement de l'eau (kg/kg)
+      REAL, INTENT(IN) :: dwatercap_dif(ngrid)  ! trend related to permanent ice
       REAL, INTENT(INOUT) :: pdqsdif(ngrid,nq)    ! tendance towards surface 
                                  !   (kg/kg.s-1)
-      REAL, INTENT(IN) :: h2oflux(ngrid)       ! value of vapour flux of H2O
-                                      ! same sign as pdqsdif
 
 c     Output:
@@ -51 +50 @@
       REAL extrasublim ! sublimation in excess of surface ice
       REAL tmpratio(ngrid)   ! D/H ratio in flux to surf
+      REAL h2oflux(ngrid)       ! value of vapour flux of H2O
+                                      ! same sign as pdqsdif
 
       INTEGER ig,l
@@ -62 +63 @@
 
         DO ig=1,ngrid
+              
+            h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice) + 
+     &          dwatercap_dif(ig)
 
             !! IF Sublimation
@@ -87 +91 @@
      &           .gt.pqsurf(ig,igcm_h2o_ice)) then
 
-C               This is how much we sublimate in excess of
+C               dwatercap_dif is how much we sublimate in excess of
 C               pqsurf for H2O                        
-                  extrasublim=
-     &              ( pqsurf(ig,igcm_h2o_ice) /ptimestep 
-     &                   +h2oflux(ig))
-c               extrasublim has same sign conv. as pdqsdif
- 
 C               hdoflux(ig) is the flux of HDO from atm. to surf.
 c               The D/H of the old ice is supposed to be 5 SMOW
 c               We need D/H of the flux to be 5, so we need
-c               extrasublim* 5 * 2 * 155.76e-6 (=1 SMOW)
+c               dwatercap_dif* 5 * 2 * 155.76e-6 (=1 SMOW)
                     hdoflux(ig)= hdoflux(ig)
-     &                       +(extrasublim*(2.*155.76e-6)*5.)
-!    &                       +extrasublim
-!                 hdoflux(ig) = h2oflux(ig) !test
-              endif
-                endif ! watercap
+     &                   +(dwatercap_dif(ig)*(2.*155.76e-6)*5.)
+                endif
+             endif ! watercap
 
             else ! condensation

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

@@ -3688 +3688 @@
            end do
      
+        call WRITEDIAGFI(ngrid,"watercap","Water ice thickness"
+     &                                  ,"kg.m-2",0,watercap)
         call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
      &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))

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

@@ -835 +835 @@
             CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
      &                      zt,zq,pqsurf,
-     &                      old_h2o_vap,pdqsdif,h2oflux,
+     &                      old_h2o_vap,pdqsdif,dwatercap_dif,
      &                      hdoflux)
             DO ig=1,ngrid
@@ -879 +879 @@
                 if(watercaptag(ig)) then
 c              dwatercap_dif same sign as pdqsdif (pos. toward ground)
-                   dwatercap_dif(ig) = -pdqsdif(ig,igcm_h2o_ice)
-     &                         - pqsurf(ig,igcm_h2o_ice)/ptimestep
+                   dwatercap_dif(ig) = pdqsdif(ig,igcm_h2o_ice)
+     &                         + pqsurf(ig,igcm_h2o_ice)/ptimestep
                    pdqsdif(ig,igcm_h2o_ice)=
      &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
@@ -892 +892 @@
                    endif
 
-                endif   
 c                 write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
                   z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
                   zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
-     $            zb(ig,2)*zc(ig,2) + (dwatercap_dif(ig)-
-c     $            zb(ig,2)*zc(ig,2) + (-
-     $            pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig)
+     $            zb(ig,2)*zc(ig,2) + 
+     $              (-pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig)
                   zq1temp(ig)=zc(ig,1)
-              endif ! if (.not.watercaptag(ig))
+                endif  !if .not.watercaptag(ig) 
+              endif ! if sublim more than surface
+
 c             Starting upward calculations for water :
               zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
@@ -925 +925 @@
 
                if(pqsurf(ig,igcm_h2o_ice)
-     &           +(dwatercap_dif(ig)+pdqsdif(ig,igcm_h2o_ice))*ptimestep
+     &           +pdqsdif(ig,igcm_h2o_ice)*ptimestep
      &           .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To
      &           pdtsrf(ig) = min(pdtsrf(ig),(To-ptsrf(ig))/ptimestep) ! ice melt case