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

Timestamp:

Nov 19, 2021, 5:00:58 PM (4 years ago)

Author:

jnaar

Message:

Correction of mass conservation bug in vdifc_mod.F : the evolution of h2o_ice_s and watercap were computed independantly in the subtimestep of surface water ice condensation / sublimation scheme. As the watercap tendency can only be negative, very specific conditions where the integrated overall tendency would be condensation but where a certain subtimestep would be subliming would result in a overall mass l
oss. As a correction, the monitoring of watercap is now done outside the subtimestep (much simpler). Local variables "zwatercap" and "zdw
atercap_dif" are now obsoletes and deleted.
JN

File:

• trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F (modified) (8 diffs)

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

@@ -131 +131 @@
 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
@@ -158 +157 @@
       REAL,INTENT(IN) :: watercap(ngrid)
       REAL,INTENT(OUT) :: dwatercap_dif(ngrid)
-      REAL :: zdwatercap_dif(ngrid)
 
 c    Subtimestep to compute h2o latent heat flux:
@@ -242 +240 @@
       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
@@ -892 +888 @@
             subtimestep = ptimestep/nsubtimestep(ig)
             ztsrf(ig)=ptsrf(ig)  !  +pdtsrf(ig)*subtimestep
-            zwatercap(ig)=watercap(ig)
 
 c           Debut du sous pas de temps
@@ -928 +923 @@
              zdqsdif(ig)=rho(ig)*dryness(ig)*zcdv(ig)
      &                      *(zq1temp(ig)-qsat(ig))
-             if ((-zdqsdif(ig)*subtimestep)
+c             write(*,*)'subliming more than available frost:  qsurf!'
+             if(.not.watercaptag(ig)) then
+               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
+c             pdqsdif > 0 : ice condensing
+c             pdqsdif < 0 : ice subliming
+c             write(*,*) "subliming more than available frost:  qsurf!"
                   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)
+     $           zb(ig,2)*zc(ig,2) +
+     $           (-zdqsdif(ig)) *subtimestep) *z1(ig)
                  zq1temp(ig)=zc(ig,1)
                endif  !if .not.watercaptag(ig) 
@@ -955 +944 @@
 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
+
+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)
+                zdtsrf(ig)=  zdqsdif(ig)*lh /pcapcal(ig)
               endif ! (latentheat_surfwater) then
 
+              DO ilay=2,nlay
+                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
+              ENDDO
 
 c             Subtimestep water budget :
@@ -971 +962 @@
               zqsurf(ig)= zqsurf(ig)+(
      &                       zdqsdif(ig))*subtimestep
-              zwatercap(ig)= zwatercap(ig)+
-     &                       zdwatercap_dif(ig)*subtimestep
 
 
@@ -981 +970 @@
      &           .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 :
-
+c             Integration of subtimestep temp and water budget :
+c             (btw could also compute the post timestep temp and ice
+c             by simply adding the subtimestep trend instead of this)
             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
+
+c if subliming more than qsurf(ice) and on watercaptag, water
+c sublimates from watercap reservoir (dwatercap_dif is <0)
+            if(watercaptag(ig)) then
+              if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
+     &           .gt.(pqsurf(ig,igcm_h2o_ice))) then
+                 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
+              endif! ((-pdqsdif(ig)*ptimestep)
+            endif !(watercaptag(ig)) then
 
            ENDDO ! of DO ig=1,ngrid