Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 2514)
+++ trunk/LMDZ.MARS/README	(revision 2515)
@@ -3366,2 +3366,13 @@
 Bugfix for streamfunction tool. Now asks the user for another file in which to
 find cv.
+
+
+== 6/05/2021 == JN
+vdifc_mod update !!
+ - vdifc now treats separately h2o_vap and other tracers
+ - water sublimation is now treated with an implicit scheme & adaptative subtimestep
+ - latent heat of sublimation is working with the implicit scheme
+ - the adaptative subtimestep is empirically calculated with a temperature variation criterion "dtmax"
+ - hdo_vap is now treated independantly after h2o_vap (should be fine !)
+
+
Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 2514)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 2515)
@@ -26,4 +26,5 @@
       use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol
       use hdo_surfex_mod, only: hdo_surfex
+c      use geometry_mod, only: longitude_deg,latitude_deg !  Joseph
       use dust_param_mod, only: doubleq, submicron, lifting
 
@@ -86,4 +87,5 @@
       integer,intent(in) :: nq 
       REAL,INTENT(IN) :: pqsurf(ngrid,nq)
+      REAL :: zqsurf(ngrid) ! temporary water tracer
       real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 
       real,intent(out) :: pdqdif(ngrid,nlay,nq) 
@@ -126,13 +128,16 @@
       REAL,SAVE :: acond,bcond
       
+c     Subtimestep & implicit treatment of water vapor
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
+      REAL zdqsdif(ngrid) ! subtimestep pdqsdif for water ice
+      REAL zwatercap(ngrid) ! subtimestep watercap for water ice
+      REAL ztsrf(ngrid) ! temporary surface temperature in tsub
+      REAL zdtsrf(ngrid) ! surface temperature tendancy in tsub
+
 c     For latent heat release from ground water ice sublimation    
-      REAL tsrf_lh(ngrid) ! temporary surface temperature
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
+      REAL tsrf_lh(ngrid) ! temporary surface temperature with lh effect
       REAL lh ! latent heat, formulation given in the Technical Document: 
               ! "Modeling water ice sublimation under Phoenix-like conditions", Montmessin et al. 2004  
-!      REAL tsrf_lw(ngrid)
-!      REAL alpha
-!      REAL T1,T2
-!      SAVE T1,T2
-!      DATA T1,T2/-604.3,1080.7/ ! zeros of latent heat equation for ice
 
 c    Tracers :
@@ -153,4 +158,11 @@
       REAL,INTENT(IN) :: watercap(ngrid)
       REAL,INTENT(OUT) :: dwatercap_dif(ngrid)
+      REAL :: zdwatercap_dif(ngrid)
+
+c    Subtimestep to compute h2o latent heat flux:
+      REAL :: dtmax = 0.5 ! subtimestep temp criterion
+      INTEGER tsub ! adaptative subtimestep (seconds)
+      REAL subtimestep !ptimestep/nsubtimestep
+      INTEGER nsubtimestep(ngrid) ! number of subtimestep (int)
 
 c    Mass-variation scheme :
@@ -226,7 +238,11 @@
       pdhdif(1:ngrid,1:nlay)=0
       pdtsrf(1:ngrid)=0
+      zdtsrf(1:ngrid)=0
       pdqdif(1:ngrid,1:nlay,1:nq)=0
       pdqsdif(1:ngrid,1:nq)=0
+      zdqsdif(1:ngrid)=0
+      zwatercap(1:ngrid)=0
       dwatercap_dif(1:ngrid)=0
+      zdwatercap_dif(1:ngrid)=0
 
 c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
@@ -788,21 +804,18 @@
 c       OU calcul de la valeur de q a la surface (water)  :
 c       ----------------------------------------
-        if (water) then 
-            call watersat(ngrid,ptsrf,pplev(1,1),qsat)
-        end if
 
 c      Inversion pour l'implicite sur q 
+c      Cas des traceurs qui ne sont pas h2o_vap
+c      h2o_vap est traite plus loin avec un sous pas de temps
+c      hdo_vap est traite ensuite car dependant de h2o_vap 
 c       --------------------------------
-        do iq=1,nq  !for all tracers including stormdust
+
+        do iq=1,nq  !for all tracers except water vapor
+          if ((.not. water).or.(.not. iq.eq.igcm_h2o_vap).or.
+     &       (.not. iq.eq.igcm_hdo_vap)) then
+
+
           zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
-
-          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:
-             zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
-             zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1)
-          else ! (re)-initialize zb(:,1)
-             zb(1:ngrid,1)=0
-          end if
+          zb(1:ngrid,1)=0
 
           DO ig=1,ngrid
@@ -822,16 +835,199 @@
           ENDDO
 
-          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
-           ! h2o vapour case (see further down).
-            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)) *z1(ig)
-            ENDDO
-         else if (hdo.and.(iq.eq.igcm_hdo_vap)) then
+            if ((iq.eq.igcm_h2o_ice) 
+     $      .or. (hdo.and.(iq.eq.igcm_hdo_ice) )) then
+
+               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)) *z1(ig)  !special case h2o_ice
+               ENDDO
+            else ! every other tracer
+               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) +
+     $             (-pdqsdif(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
+               ENDDO
+            endif !((iq.eq.igcm_h2o_ice)
+c         Starting upward calculations for simple mixing of tracer (dust)
+          DO ig=1,ngrid
+             zq(ig,1,iq)=zc(ig,1)
+             DO ilay=2,nlay
+               zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
+             ENDDO
+          ENDDO
+          endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then
+        enddo ! of do iq=1,nq
+
+c --------- h2o_vap --------------------------------
+
+
+c      Traitement de la vapeur d'eau h2o_vap
+c      Utilisation d'un sous pas de temps afin
+c      de decrire le flux de chaleur latente 
+
+
+        do iq=1,nq  
+          if ((water).and.(iq.eq.igcm_h2o_vap)) then
+
+
+           DO ig=1,ngrid
+             zqsurf(ig)=pqsurf(ig,igcm_h2o_ice)
+           ENDDO ! ig=1,ngrid
+
+c          make_tsub : sous pas de temps adaptatif
+c          la subroutine est a la fin du fichier
+
+           call make_tsub(ngrid,pdtsrf,zqsurf,
+     &                    ptimestep,dtmax,watercaptag,
+     &                    nsubtimestep)
+
+c           Calculation for turbulent exchange with the surface (for ice)
+c           initialization of ztsrf, which is surface temperature in
+c           the subtimestep.
+           DO ig=1,ngrid
+            subtimestep = ptimestep/nsubtimestep(ig)
+            ztsrf(ig)=ptsrf(ig)  !  +pdtsrf(ig)*subtimestep
+            zwatercap(ig)=watercap(ig)
+
+c           Debut du sous pas de temps
+
+            DO tsub=1,nsubtimestep(ig)
+
+c           C'est parti ! 
+
+             zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
+     &                     /float(nsubtimestep(ig))
+             zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
+     &                     /float(nsubtimestep(ig))
+             zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1)
+             
+             z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+             zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
+             zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+             DO ilay=nlay-1,2,-1
+                 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+                 zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
+     $            zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+                 zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+             ENDDO  
+             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)) * z1(ig)
+
+             call watersat(1,ztsrf(ig),pplev(1,1),qsat(ig))
+             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)
+
+             zdqsdif(ig)=rho(ig)*dryness(ig)*zcdv(ig)
+     &                      *(zq1temp(ig)-qsat(ig))
+             if ((-zdqsdif(ig)*subtimestep)
+     &            .gt.(zqsurf(ig))) then
+c             write(*,*)'subliming more than available frost:  qsurf!'
+               if(watercaptag(ig)) then
+c             dwatercap_dif same sign as pdqsdif (pos. toward ground)
+                  zdwatercap_dif(ig) = zdqsdif(ig)
+     &                        + zqsurf(ig)/subtimestep
+                  zdqsdif(ig)=
+     &                        -zqsurf(ig)/subtimestep
+c                  write(*,*) "subliming more than available frost:  qsurf!"
+               else  !  (case not.watercaptag(ig)) 
+                  zdqsdif(ig)=
+     &                        -zqsurf(ig)/subtimestep
+                  zdwatercap_dif(ig) = 0.0
+
+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) + 
+     $           (-zdqsdif(ig)-zdwatercap_dif(ig)) *subtimestep) *z1(ig)
+                 zq1temp(ig)=zc(ig,1)
+               endif  !if .not.watercaptag(ig) 
+             endif ! if sublim more than surface
+
+c             Starting upward calculations for water :
+c             Actualisation de h2o_vap dans le premier niveau
+             zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
+              
+c             Take into account the H2O latent heat impact on the surface temperature
+              if (latentheat_surfwater) then
+                lh=(2834.3-0.28*(ztsrf(ig)-To)-
+     &              0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3
+                zdtsrf(ig)=  (zdwatercap_dif(ig)+
+     &               zdqsdif(ig))*lh /pcapcal(ig)
+              endif ! (latentheat_surfwater) then
+
+
+c             Subtimestep water budget :
+
+              ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig) + zdtsrf(ig))
+     &                          *subtimestep
+              zqsurf(ig)= zqsurf(ig)+(
+     &                       zdqsdif(ig))*subtimestep
+              zwatercap(ig)= zwatercap(ig)+
+     &                       zdwatercap_dif(ig)*subtimestep
+
+
+c             We ensure that surface temperature can't rise above the solid domain if there
+c             is still ice on the surface (oldschool)
+               if(zqsurf(ig)
+     &           +zdqsdif(ig)*subtimestep
+     &           .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To
+     &           zdtsrf(ig) = min(zdtsrf(ig),(To-ztsrf(ig))/subtimestep) ! ice melt case
+     
+
+              DO ilay=2,nlay
+                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
+              ENDDO
+
+c             Fin du sous pas de temps
+
+            ENDDO ! tsub=1,nsubtimestep 
+
+c             Integration of subtimestep water budget :
+
+            pdtsrf(ig)= (ztsrf(ig) -
+     &                     ptsrf(ig))/ptimestep
+            pdqsdif(ig,igcm_h2o_ice)=
+     &                  (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice))/ptimestep
+            dwatercap_dif(ig)=(zwatercap(ig) -
+     &                     watercap(ig))/ptimestep
+
+           ENDDO ! of DO ig=1,ngrid
+          END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
+
+c --------- end of h2o_vap ----------------------------
+
+c --------- hdo_vap -----------------------------------
+
+c    hdo_ice has already been with along h2o_ice
+c    amongst "normal" tracers (ie not h2o_vap)
+
+         if (hdo.and.(iq.eq.igcm_hdo_vap)) then
+          zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
+          zb(1:ngrid,1)=0
+
+          DO ig=1,ngrid
+               z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
+               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+          ENDDO
+
+          DO ilay=nlay-1,2,-1
+               DO ig=1,ngrid
+                z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $           zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+                zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
+     $           zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+                zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+               ENDDO
+          ENDDO
+
             CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
      &                      zt,zq,pqsurf,
@@ -846,103 +1042,19 @@
             ENDDO
 
-          else ! general case
             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) +
-     $        (-pdqsdif(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
+              zq(ig,1,iq)=zc(ig,1)
+              DO ilay=2,nlay
+                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
+              ENDDO
             ENDDO
-
-          endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
-  
-          IF ((water).and.(iq.eq.igcm_h2o_vap)) then 
-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)
-
-              pdqsdif(ig,igcm_h2o_ice)=rho(ig)*dryness(ig)*zcdv(ig)
-     &                       *(zq1temp(ig)-qsat(ig))
-c             write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
-            END DO
-
-            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
-c              write(*,*)'subliming more than available frost:  qsurf!'
-                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
-                   pdqsdif(ig,igcm_h2o_ice)=
-     &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
-c                  write(*,*) "subliming more than available frost:  qsurf!"
-                else  !  (case not.watercaptag(ig)) 
-                   pdqsdif(ig,igcm_h2o_ice)=
-     &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
-                   dwatercap_dif(ig) = 0.0
-                   if (hdo) then
-                     h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice)
-                   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) + 
-     $              (-pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig)
-                  zq1temp(ig)=zc(ig,1)
-                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)
-              
-!c             Take into account H2O latent heat in surface energy budget
-!c             We solve dT/dt = (2834.3-0.28*(T-To)-0.004*(T-To)^2)*1e3*iceflux/cpp
-!              tsrf_lw(ig) = ptsrf(ig) + pdtsrf(ig) *ptimestep
-!             
-!              alpha = exp(-4*abs(T1-T2)*pdqsdif(ig,igcm_h2o_ice)
-!     &            *ptimestep/pcapcal(ig))
-!
-!              tsrf_lw(ig) = (tsrf_lw(ig)*(T2-alpha*T1)+T1*T2*(alpha-1))
-!     &                      /(tsrf_lw(ig)*(1-alpha)+alpha*T2-T1)  ! surface temperature at t+1
-!
-!              pdtsrf(ig) = (tsrf_lw(ig)-ptsrf(ig))/ptimestep
-c             Take into account the H2O latent heat impact on the surface temperature
-              if (latentheat_surfwater) then
-                tsrf_lh(ig) = ptsrf(ig) + pdtsrf(ig) *ptimestep
-                lh=(2834.3-0.28*(tsrf_lh(ig)-To)-
-     &              0.004*(tsrf_lh(ig)-To)*(tsrf_lh(ig)-To))*1.e+3
-                pdtsrf(ig)= pdtsrf(ig) + (dwatercap_dif(ig)+
-     &               pdqsdif(ig,igcm_h2o_ice))*lh /pcapcal(ig)
-              endif
-
-               if(pqsurf(ig,igcm_h2o_ice)
-     &           +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
-     
-            ENDDO ! of DO ig=1,ngrid
-          ELSE
-c           Starting upward calculations for simple mixing of tracer (dust)
-            DO ig=1,ngrid
-               zq(ig,1,iq)=zc(ig,1)
-            ENDDO
-          END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
-
-          DO ilay=2,nlay
-             DO ig=1,ngrid
-               zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
-             ENDDO
-          ENDDO
+         endif ! (hdo.and.(iq.eq.igcm_hdo_vap))
+
+c --------- end of hdo ----------------------------
+
         enddo ! of do iq=1,nq
       end if ! of if(tracer)
+
+c --------- end of tracers  ----------------------------
+
 
 C       Diagnostic output for HDO
@@ -1004,4 +1116,35 @@
 c====================================
 
-
+      SUBROUTINE make_tsub(naersize,dtsurf,qsurf,ptimestep,
+     $                     dtmax,watercaptag,ntsub)
+
+c Pas de temps adaptatif en estimant le taux de sublimation
+c et en adaptant avec un critere "dtmax" du chauffage a accomoder
+c dtmax est regle empiriquement (pour l'instant) a 0.5 K
+
+      integer,intent(in) :: naersize
+      real,intent(in) :: dtsurf(naersize)
+      real,intent(in) :: qsurf(naersize)
+      logical,intent(in) :: watercaptag(naersize)
+      real,intent(in) :: ptimestep
+      real,intent(in)  :: dtmax
+      real  :: ztsub(naersize)
+      integer  :: i
+      integer,intent(out) :: ntsub(naersize)
+
+      do i=1,naersize
+        if ((qsurf(i).eq.0).and.
+     &      (.not.watercaptag(i))) then
+          ntsub(i) = 1
+        else 
+          ztsub(i) = ptimestep * dtsurf(i) / dtmax
+          ntsub(i) = ceiling(abs(ztsub(i)))
+        endif ! (qsurf(i).eq.0) then
+c      
+c       write(78,*), dtsurf*ptimestep, dtsurf, ntsub
+      enddo! 1=1,ngrid
+
+
+
+      END SUBROUTINE make_tsub
       END MODULE vdifc_mod
