Changeset 2947 for trunk/LMDZ.MARS/libf/phymars/dyn1d

Timestamp:

Apr 26, 2023, 11:29:35 AM (20 months ago)

Author:

llange

Message:

Mars PCM
Following r-2942: Fix a bug in the 1D: inertiesoil was not written in
the startfile
I introduce the possibility to prescribe subsurface ice (when no
startfi file): In the callphys.def, set subsurface_ice_depth to the
depth where you want ice (set value <0 for no ice). It will assign for
these depth a thermal inertia of 2100.
LL

File:

• trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F (modified) (4 diffs)

trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F

@@ -158 +158 @@
       LOGICAL :: found
       REAL albedo_read(1,2,1)      ! surface albedo
+
+c   LL: Possibility to add subsurface ice
+      REAL :: ice_depth            ! depth of the ice table, ice_depth < 0. means no ice
+      REAL :: inertieice = 2100.   ! ice thermal inertia
+      integer :: iref
+ 
 c=======================================================================
 
@@ -528 +534 @@
       write(*,*) " inertia = ",inertiedat(1,1)
 
+      ice_depth = -1 ! default value: no ice
+      CALL getin("subsurface_ice_depth",ice_depth)
+
       z0(1)=z0_default ! default value for roughness
       write(*,*) 'Surface roughness length z0 (m)?'
@@ -690 +699 @@
 ! ------------------------------------------
       volcapa=1.e6 ! volumetric heat capacity
+
       if(.not. startfile_1D ) then
-      DO isoil=1,nsoil
-         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
+
+! Initialize depths
+! -----------------
+       do isoil=1,nsoil
+         layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
+       enddo
+
+! Creating the new soil inertia table if there is subsurface ice :
+       IF (ice_depth.gt.0) THEN 
+         iref = 1 ! ice/regolith boundary index
+           IF (ice_depth.lt.layer(1)) THEN
+             inertiedat(1,1) = sqrt( layer(1) /
+     &        ( (ice_depth/inertiedat(1,1)**2) +
+     &        ((layer(1)-ice_depth)/inertieice**2) ) )
+             DO isoil=1,nsoil 
+              inertiedat(1,isoil) = inertieice
+             ENDDO    
+           ELSE ! searching for the ice/regolith boundary :
+           DO isoil=1,nsoil 
+            IF ((ice_depth.ge.layer(isoil)).and.
+     &         (ice_depth.lt.layer(isoil+1))) THEN
+                  iref=isoil+1
+                  EXIT
+            ENDIF
+           ENDDO
+!         We then change the soil inertia table :
+           DO isoil=1,iref-1 
+            inertiedat(1,isoil) = inertiedat(1,1)
+           ENDDO
+!         We compute the transition in layer(iref)
+           inertiedat(1,iref) = sqrt( (layer(iref)-layer(iref-1)) /
+     &          ( ((ice_depth-layer(iref-1))/inertiedat(1,1)**2) +
+     &          ((layer(iref)-ice_depth)/inertieice**2) ) ) 
+!         Finally, we compute the underlying ice :
+           DO isoil=iref+1,nsoil
+            inertiedat(1,isoil) = inertieice
+           ENDDO
+         ENDIF ! (ice_depth.lt.layer(1))      
+       ELSE ! ice_depth <0 all is set to surface thermal inertia
+         DO isoil=1,nsoil
+           inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
+         ENDDO
+       ENDIF ! ice_depth.gt.0
+
+       inertiesoil(1,:,1) = inertiedat(1,:)
+
+       DO isoil = 1,nsoil
          tsoil(isoil)=tsurf(1)  ! soil temperature
-      ENDDO
+       ENDDO
+
       endif !(.not. startfile_1D )
 
@@ -789 +845 @@
       call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
      &              dtphys,time,
-     &              tsurf,tsoil,albedo,emis,q2,qsurf,tauscaling,
+     &              tsurf,tsoil,inertiesoil,albedo,emis,
+     &              q2,qsurf,tauscaling,
      &              totcloudfrac,wstar,watercap)
       endif !(.not. startfile_1D )