Changeset 4119 for trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90

Timestamp:

Mar 10, 2026, 10:58:47 PM (5 weeks ago)

Author:

tbertrand

Message:

Pluto PCM:
Some fixing for the paleo mode
TB

File:

• trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90 (modified) (13 diffs)

trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90

@@ -48 +48 @@
       REAL tgla(ngrid),tbase(ngrid)   !temperature at the base of glacier different of ptsurf
       REAL :: zdqsurf(ngrid)   ! tendency of qsurf
+      REAL :: diff(ngrid)   ! height difference (topography+ice layer) between two adjacent points
       REAL massflow  ! function
-      REAL hlim,qlim,difflim,diff,stock,H0,totmasstmp
-      INTEGER compt !compteur
+      REAL hlim,qlim,difflim,stock,H0,totmasstmp
+      INTEGER compt,nbedge 
       INTEGER slid ! option slid 0 or 1
       INTEGER :: edge(4) ! index of the adjecent points, N,S,E,W
@@ -63 +64 @@
       difflim=0.5 ! limit height (m) between two adjacent point to start spreading the glacier
       zdqsurf(:)=0.
+
+      ! We set the number of edges to make it applicable for 1D case
+      IF (iim.eq.1) THEN
+         nbedge=2   ! 1D Case: Only check North/South
+      ELSE
+         nbedge=4   ! 3D Cse: Check N,S,E,W
+      ENDIF
+
 !--------------- Dimensions -------------------------------------
       
@@ -73 +82 @@
       qlim=hlim*1000. ! kg
       !! Security for not depleted all ice too fast in one timestep 
-      secu=4 
-
+      secu=2 
+ 
       !*************************************
       ! Loop over all points 
       !*************************************
       DO ig=1,ngrid
-         !if (ig.eq.ngrid) then
-         !  print*, 'qpole=',pqsurf(ig,igcm_n2),qlim
-         !endif
 
          !*************************************
          ! if significant amount of ice, where qsurf > qlim 
          !*************************************
+
          IF (pqsurf(ig,igcm_n2).gt.qlim) THEN
 
@@ -94 +101 @@
               tgla(ig)=ptsurf(ig)
               ! temperature at the base (we neglect convection)
-              tbase(ig)=tgla(ig)+20.*pqsurf(ig,igcm_n2)/1.e6 ! Umurhan et al.
-              if (tbase(ig).gt.63.) then  
+              ! tbase(ig)=tgla(ig)+20.*pqsurf(ig,igcm_n2)/1.e6 ! Umurhan et al.
+
+              ! DEBUG test: paleo simulation without basal melting
+              tbase(ig) = tgla(ig)
+
+              if (tbase(ig).gt.63.) then 
                  slid=1   ! activate slide
               else
@@ -150 +161 @@
 
               masstmp(:)=0. ! mass to be transferred
+              diff(:)=0. ! height difference between adjacent points (m)
               totmasstmp=0. ! total mass to be transferred
               H0=phisfi0(ig)/g+pqsurf(ig,igcm_n2)/1000. ! height of ice at ig (m)
@@ -161 +173 @@
                  !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point
                  DO i=inddeb,indfin
-                    diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.)  ! Height difference between ig and adjacent points (m)
-                    if (diff.gt.difflim) then
-                     masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s)
+                    diff(i)=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.)  ! Height difference between ig and adjacent points (m)
+                    if (diff(i).gt.difflim) then
+                     masstmp(i)=massflow(ig,i,tgla(ig),diff(i),pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s)
                     else
                      masstmp(i)=0.
@@ -171 +183 @@
                  if (totmasstmp.gt.0.) then
                    !! 2) Available mass to be transferred
-                   stock=maxval(masstmp(:))/secu ! kg/m2/s
+                   stock=min(maxval(diff(:)*1000./timstep),maxval(masstmp(:))/secu) ! kg/m2/s
                    !! 3) Real amounts of ice to be transferred :
                    masstmp(:)=masstmp(:)/totmasstmp*stock*cell_area(ig)/cell_area(inddeb)  !kg/m2/s   all area around the pole are the same
@@ -217 +229 @@
 
                  !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point
-                 DO compt=1,4
+                 DO compt=1,nbedge
                     i=edge(compt)
-                    diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.)  ! (m)
-                    if (diff.gt.difflim) then
-                     masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid)
+                    diff(i)=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.)  ! (m)
+
+
+                    if (diff(i).gt.difflim) then
+                     masstmp(i)=massflow(ig,i,tgla(ig),diff(i),pqsurf(ig,igcm_n2),timstep,slid)
+                    
+
                     else
                      masstmp(i)=0.
@@ -229 +245 @@
                  if (totmasstmp.gt.0) then
                    !! 2) Available mass to be transferred
-                   stock=maxval(masstmp(:))/secu ! kg/m2/s
+                   stock=min(maxval(diff(:)*1000./timstep),maxval(masstmp(:))/secu) ! kg/m2/s
                    !! 3) Real amounts of ice to be transferred :
-                   DO compt=1,4
+                   DO compt=1,nbedge
                     i=edge(compt)
                     masstmp(i)=masstmp(i)/totmasstmp*stock*cell_area(ig)/cell_area(i)  !kg/m2/s
@@ -246 +262 @@
                      pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4
                      !! add CH4
-                     DO compt=1,4
+                     DO compt=1,nbedge
                        i=edge(compt)
                        pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 
@@ -257 +273 @@
                      pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco
                      !! add CO
-                     DO compt=1,4
+                     DO compt=1,nbedge 
                        i=edge(compt)
                        pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco 
@@ -265 +281 @@
                    ! Add N2
                    totmasstmp=0.
-                   DO compt=1,4
+                   DO compt=1,nbedge
                        i=edge(compt)
                        pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep
@@ -285 +301 @@
                 
       ENDDO  ! ig 
+
 
       end subroutine spreadglacier_paleo