Changeset 3098

Timestamp:

Oct 23, 2023, 7:02:06 PM (14 months ago)

Author:

evos

Message:

Mars PCM

We added the to module vdifc the possibily of subsurface intercation, mostly to have the option of buried glaicer that can lose ice and create polar layers.
if no ice is on the surface the ice the atmosphere can directly interact with the subsurface ice.
if ice is on the surface it can interact with the subsurface ice.
For now it is not working with more then one subslope. I will add this feature soon
In order to use it there are two flags that needs to be set to true.
paleoclimate=.true.
lag_layer=.true.
in addition one need to give the following parametrs in startfi or run.def:
d_coef=X for the diffusion coeficent (4e-4 is the defualt)
<h2o_ice_depth=X> if its zero or lower there is no subsurface ice (default -1).
outputs: zdq_ssi, momentarily zdqsdif after the interaction with the SSIzdq_subtimestep,  momentarily flux after the interaction with the SSIzdq_ssi_frost momentarily zdqsdif after the interaction of the frost with the SSI
q1, atmospheric first layer water vapor quatity
qeq, SSI water vapor quantity

Eran Vos

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/conf_phys.F (modified) (2 diffs)
• libf/phymars/dyn1d/init_testphys1d_mod.F90 (modified) (2 diffs)
• libf/phymars/dyn1d/testphys1d.F90 (modified) (1 diff)
• libf/phymars/paleoclimate_mod.F90 (modified) (3 diffs)
• libf/phymars/phyetat0_mod.F90 (modified) (4 diffs)
• libf/phymars/phyredem.F90 (modified) (2 diffs)
• libf/phymars/physiq_mod.F (modified) (1 diff)
• libf/phymars/vdifc_mod.F (modified) (20 diffs)

trunk/LMDZ.MARS/changelog.txt

@@ -4273 +4273 @@
 Prettyfied solarlong.F and made it a module. Likewise for conf_phys.F
 
+== 23/10/2023 == EV
+We added the to module vdifc the possibily of subsurface intercation, mostly to have the option of buried glaicer that can lose ice and create polar layers.
+if no ice is on the surface the ice the atmosphere can directly interact with the subsurface ice.
+if ice is on the surface it can interact with the subsurface ice.
+**For now it is not working with more then one subslope. I will add this feature soon**
+In order to use it there are two flags that needs to be set to true.
+paleoclimate=.true.
+lag_layer=.true.
+in addition one need to give the following parametrs in startfi or run.def:
+d_coef=X for the diffusion coeficent (4e-4 is the defualt)
+<h2o_ice_depth=X> if its zero or lower there is no subsurface ice (default -1).
+outputs: zdq_ssi, momentarily zdqsdif after the interaction with the SSIzdq_subtimestep,  momentarily flux after the interaction with the SSIzdq_ssi_frost momentarily zdqsdif after the interaction of the frost with the SSI
+q1, atmospheric first layer water vapor quatity
+qeq, SSI water vapor quantity
+
+
 Some code tidying:
 Made pi in module comcstfi_h a parameter (and not a variable recomputed at

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

@@ -37 +37 @@
       use aeropacity_mod, only: iddist, topdustref
       USE mod_phys_lmdz_transfert_para, ONLY: bcast
-      USE paleoclimate_mod,ONLY: paleoclimate,albedo_perenialco2
+      USE paleoclimate_mod,ONLY: paleoclimate,albedo_perenialco2,
+     &                           lag_layer 
       use microphys_h, only: mteta
 
@@ -338 +339 @@
 ! PALEOCLIMATE
 
-         write(*,*)"Is it paleoclimate run?"
+         write(*,*)"Using lag layer??"
+         lag_layer=.false.
+         call getin_p("lag_layer",lag_layer)
+         write(*,*) " lag_layer = ", lag_layer
+
+        write(*,*)"Is it paleoclimate run?"
          paleoclimate=.false. ! default value
          call getin_p("paleoclimate",paleoclimate)
          write(*,*)" paleoclimate = ",paleoclimate
+
 
          write(*,*)"Albedo for perenial CO2 ice?"

trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

@@ -477 +477 @@
     call getin("subsurface_ice_depth",ice_depth)
 
+
     z0(1) = z0_default ! default value for roughness
     write(*,*) 'Surface roughness length z0 (m)?'
@@ -716 +717 @@
         write(*,*) 'Prescribed atmospheric water vapor profile'
         write(*,*) 'Unless it reaches saturation (maximal value)'
+    else if (atm_wat_profile .eq. 2) then
+        write(*,*) 'Prescribed atmospheric water vapor profile'
+        write(*,*) 'like present day northern midlatitude'
+        write(*,*) 'Unless it reaches saturation (maximal value)'
     else
         error stop 'Water vapor profile value not correct!'

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

@@ -179 +179 @@
             endif
         endif
-    endif
+    !EV profile
+!            IF(atm_wat_profile.eq.2) THEN
+!             DO ilayer=1,16
+!              q(1,ilayer,igcm_h2o_vap)=min(zqsat(ilayer),6e-5)
+!             ENDDO! ilayer=1,16
+!             DO ilayer=17,22
+!              q(1,ilayer,igcm_h2o_vap)=min(zqsat(ilayer),3e-5)
+!             ENDDO! ilayer=17,22
+!             DO ilayer=23,nlayer
+!              q(1,ilayer,igcm_h2o_vap)=0
+!             ENDDO! ilayer=23,nlayer
+!             endif
+     endif
 
     ! Call physics

trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90

@@ -15 +15 @@
 !$OMP THREADPRIVATE(paleoclimate)
 
-    real, save, allocatable :: lag_h2o_ice(:,:)  ! Thickness of the lag before H2O ice [m]
+    real, save, allocatable :: h2o_ice_depth(:,:)  ! Thickness of the lag before H2O ice [m]
     real, save, allocatable :: lag_co2_ice(:,:)  ! Thickness of the lag before CO2 ice [m]
     real, save :: albedo_perenialco2             ! Albedo for perenial co2 ice [1]
-!$OMP THREADPRIVATE(lag_h2o_ice,lag_co2_ice,albedo_perenialco2)
+    real, save, allocatable :: d_coef(:,:)  ! Diffusion coeficent
+    LOGICAL,SAVE :: lag_layer ! does lag layer is present?
+!$OMP THREADPRIVATE(h2o_ice_depth,lag_co2_ice,albedo_perenialco2)
 
     CONTAINS
@@ -29 +31 @@
   integer,intent(in) :: nslope ! number of slope within a mesh 
 
-    allocate(lag_h2o_ice(ngrid,nslope))
+    allocate(h2o_ice_depth(ngrid,nslope))
     allocate(lag_co2_ice(ngrid,nslope))
+    allocate(d_coef(ngrid,nslope))
   end subroutine ini_paleoclimate_h
 
@@ -36 +39 @@
 
   implicit none
-    if (allocated(lag_h2o_ice)) deallocate(lag_h2o_ice)
+    if (allocated(d_coef)) deallocate(d_coef)
+    if (allocated(h2o_ice_depth)) deallocate(h2o_ice_depth)
     if (allocated(lag_co2_ice)) deallocate(lag_co2_ice)
   end subroutine end_paleoclimate_h

trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90

@@ -27 +27 @@
   use comsoil_h, only: flux_geo
   USE comslope_mod, ONLY: nslope, major_slope
-  USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice
+  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef
   USE comcstfi_h, only: pi
   use geometry_mod, only: latitude
@@ -744 +744 @@
   if (startphy_file) then
 ! Depth of H2O lag
-   call get_field("lag_h2o_ice",lag_h2o_ice,found,indextime)
-   if (.not.found) then
-     write(*,*) "phyetat0: Failed loading <lag_h2o_ice> : ", &
-                          "<lag_h2o_ice> is set as -1 (no subsurface ice)"
-     lag_h2o_ice(:,:) = -1.
-   endif
+   call get_field("h2o_ice_depth",h2o_ice_depth,found,indextime)
+   if (.not.found) then
+     write(*,*) "phyetat0: Failed loading <h2o_ice_depth> : ", &
+                          "<h2o_ice_depth> is set as -1 (no subsurface ice)"
+     h2o_ice_depth(:,:) = -1.
+   endif
+
+! Diffuesion coeficent 
+   call get_field("d_coef",d_coef,found,indextime)
+   if (.not.found) then
+     write(*,*) "phyetat0: Failed loading <d_coef> : ", &
+                          "<d_coef> is set as 4e-4 (defualt value)"
+     d_coef(:,:) = 4e-4
+   endif
+
 
 !  Depth of CO2 lag
@@ -774 +783 @@
    endif ! notfound
   else ! no startfiphyle
-     lag_h2o_ice(:,:) = -1.
+     h2o_ice_depth(:,:) = -1.
      lag_co2_ice(:,:) = -1.
      perenial_co2ice(:,:) = 0.
+     d_coef(:,:)= 4.e-4
      if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
        DO islope = 1,nslope
@@ -786 +796 @@
   endif !startphy_file
 else
-  lag_h2o_ice(:,:) = -1.
+  h2o_ice_depth(:,:) = -1.
   lag_co2_ice(:,:) = -1.
   perenial_co2ice(:,:) = 0.
+  d_coef(:,:)= 4.e-4
 endif !paleoclimate
 

trunk/LMDZ.MARS/libf/phymars/phyredem.F90

@@ -25 +25 @@
   use time_phylmdz_mod, only: daysec
   use comslope_mod, ONLY: nslope
-  USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice  
+  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice  
   implicit none
  
@@ -156 +156 @@
   ! Paleoclimate outputs
   if(paleoclimate) then
-    call put_field("lag_h2o_ice","Depth of the H2O lags",lag_h2o_ice)
+    call put_field("h2o_ice_depth","Depth to the fisrt H2O ice",h2o_ice_depth)
     call put_field("lag_co2_ice","Depth of the CO2 lags",lag_co2_ice)
   endif

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

@@ -1514 +1514 @@
      &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
      &        zcondicea_co2microp,sensibFlux,
-     &        dustliftday,local_time,watercap,dwatercap_dif)
+     &        dustliftday,local_time,watercap,dwatercap_dif,
+     &        tsoil,nsoilmx)
 
           DO ig=1,ngrid

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

@@ -13 +13 @@
      $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
      $                hfmax,pcondicea_co2microp,sensibFlux,
-     $                dustliftday,local_time,watercap, dwatercap_dif)
-
+     $                dustliftday,local_time,watercap, dwatercap_dif,
+     $                tsoil,nsoil)
       use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
      &                      igcm_dust_submicron, igcm_h2o_vap,
@@ -32 +32 @@
      &                      subslope_dist,major_slope,iflat
       use microphys_h, only: To
+      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
+      use comsoil_h, only: layer, mlayer
 
       IMPLICIT NONE
@@ -61 +63 @@
 c   ----------
 
-      INTEGER,INTENT(IN) :: ngrid,nlay
+      INTEGER,INTENT(IN) :: ngrid,nlay,nsoil
       REAL,INTENT(IN) :: ptimestep
       REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
@@ -75 +77 @@
       REAL,INTENT(IN) :: pcapcal(ngrid,nslope)
       REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
+      REAL,INTENT(IN) :: tsoil(ngrid,nsoil,nslope)
 
 c    Argument added for condensation:
@@ -102 +105 @@
       REAL :: pt(ngrid,nlay)
  
-      INTEGER ilev,ig,ilay,nlev,islope
+      INTEGER ilev,ig,ilay,nlev,islope,ik,lice
 
       REAL z4st,zdplanck(ngrid)
@@ -118 +121 @@
       REAL zcst1
       REAL zu2(ngrid)
+      REAL Tice(ngrid,nslope) ! subsurface temperature where ice is located.
+      REAL qeq(ngrid,nslope) ! saturation water vapor in the subsurface
+      REAL dist_up(ngrid,nslope) !distance from ice to layer above
+      REAL dist_down(ngrid,nslope) !distance from ice to layer down
+      REAL dist_sum(ngrid,nslope) ! sum of distance
+      REAL zdqsdif_ssi(ngrid,nslope) !SSI flux
+      REAL zdqsdif_ssi_frost(ngrid,nslope) !SSI-frost interaction
 
       EXTERNAL SSUM,SCOPY
@@ -155 +165 @@
       REAL rho(ngrid) ! near surface air density
       REAL qsat(ngrid)
+      REAL qsat2(ngrid,nslope)
+      REAL resist(ngrid,nslope) !subsurface ice flux reduction coef
 
       REAL hdoflux(ngrid,nslope)       ! value of vapour flux of HDO
@@ -272 +284 @@
       zdqsdif(1:ngrid)=0
       dwatercap_dif(1:ngrid,1:nslope)=0
-
+!      h2o_ice_depth(1:ngrid,1:nslope)=5
 c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
 c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
@@ -914 +926 @@
 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
@@ -920 +932 @@
           DO islope = 1,nslope
            DO ig=1,ngrid
+
+             
+
              zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/
      &             cos(pi*def_slope_mean(islope)/180.)
@@ -943 +958 @@
             zq_tmp_vap(ig,:,:) =zq(ig,:,:)
 c           Debut du sous pas de temps
-
             DO tsub=1,nsubtimestep(ig)
 
@@ -980 +994 @@
                if ((-zdqsdif(ig)*subtimestep)
      &            .gt.(zqsurf(ig))) then
+                
+    
+                 !EV subsurface ice
+                IF(h2o_ice_depth(ig,islope) .gt. 0 .and. lag_layer)
+     &           then 
+                zdqsdif(ig)=
+     &                        -zqsurf(ig)/subtimestep
+                zqsurf(ig)=0
+
+               DO ik=0,nsoil-2 ! go through all the layers to find the ice locations
+                IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.
+     &           (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN
+                 lice = ik+1
+                EXIT
+                ENDIF
+               ENDDO !of subsurface loop
+                IF (lice .gt. 1) then !calculate the distance from the layers
+                 dist_up(ig,islope)=(h2o_ice_depth(ig,islope)
+     &                              -mlayer(lice-1))
+                 dist_down(ig,islope)=(mlayer(lice)
+     &                                -h2o_ice_depth(ig,islope))
+                 dist_sum(ig,islope)=dist_up(ig,islope)
+     &                               +dist_down(ig,islope)
+                 Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
+     &                          *tsoil(ig,lice-1,islope)
+     &           /dist_sum(ig,islope))+
+     &           (dist_down(ig,islope)*tsoil(ig,lice,islope)
+     &           /dist_sum(ig,islope))
+                        ELSE
+                                Tice(ig,islope)=tsoil(ig,1,islope)
+                        ENDIF
+              call watersat(1,Tice(ig,1),pplev(ig,1)
+     &                     ,qsat2(ig,1))
+              qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))
+     &                        *qsat2(ig,1)
+              resist(ig,1)=(1+(h2o_ice_depth(ig,islope)*zcdh(ig)
+     &                           /d_coef(ig,islope)))
+              !write(*,*)'R=',resist(ig,islope)
+              !write(*,*)'zice=',h2o_ice_depth(ig,islope)
+              !write(*,*)'D=',d_coef(ig,islope)
+              !write(*,*)'zcdh=',zcdh(ig)
+              zb(ig,1)=zb(ig,1)/resist(ig,1) ! change zb to account subsurface ice
+              !vdifc algorithem !!!!!needs to change  reseist io to the mean
+              !beacuse of the slopes!!!!!
+              z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+             zc(ig,nlay)=za(ig,nlay)*zq_tmp_vap(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_tmp_vap(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_tmp_vap(ig,1,iq)+
+     $        zb(ig,2)*zc(ig,2)) * z1(ig)
+             zd(ig,1)=zb(ig,1)*z1(ig)
+             zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
+             zdqsdif_ssi(ig,1)=rho(ig)*dryness(ig)*zcdv(ig)    
+     &                      *(zq1temp(ig)-qeq(ig,islope))
+                      call write_output('zdq_ssi',
+     &                '','',zdqsdif_ssi(ig,1))
+                       call write_output('qeq',
+     &                '','',qeq(ig,1))
+                       call write_output('q1',
+     &                '','',zq1temp(ig))
+             !write(*,*)'qeq=',qeq(ig,1)
+             !write(*,*)'q1=',zq1temp(ig) 
+             !write(*,*)'zdqsdif_ssi=',zdqsdif_ssi(ig,1)
+             !I should some all the interactions with the SSI accoridng
+             !to the statistics and the evaluate
+             
+             !write(*,*) "zdq_ssi*t", zdqsdif_ssi(ig)*subtimestep
+             if (zdqsdif_ssi(ig,1)<0) then
+               !zdqsdif(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep)
+               zdqsdif(ig)=zdqsdif(ig)+zdqsdif_ssi(ig,1)
+                call write_output('zdq_zdqssi',
+     &                '','',zdqsdif(ig)+zdqsdif_ssi(ig,1))
+             endif
+              ELSE
 c             pdqsdif > 0 : ice condensing
 c             pdqsdif < 0 : ice subliming
-c             write(*,*) "subliming more than available frost:  qsurf!"
+             write(*,*) "subliming more than available frost:  qsurf!"
                   zdqsdif(ig)=
      &                        -zqsurf(ig)/subtimestep
@@ -991 +1087 @@
      $           (-zdqsdif(ig)) *subtimestep) *z1(ig)
                  zq1temp(ig)=zc(ig,1)
+                 ENDIF !if h2o_ice_depth>0 and lag_layer
                endif  !if .not.watercaptag(ig) 
              endif ! if sublim more than surface
+
+              ! subsurface ice <--> frost interaction !EV
+             IF(h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
+     &        .and. zqsurf(ig) .gt. 0) then
+               zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
+     &                         /h2o_ice_depth(ig,1))
+     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
+               !needs to change to the mean of eq
+               zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1) 
+                !!!!! zdsqdif_ssi_frosst need to be changed to an
+                !average
+              ELSEIF (h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
+     &        .and. watercaptag(ig)) then
+              zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
+     &                         /h2o_ice_depth(ig,1))
+     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
+               zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1)
+               !needs to change to the mean of eq
+              ENDIF
+              call write_output('zdqsdif_ssi_frost',
+     &                '','',zdqsdif_ssi_frost(ig,1))
+              call write_output('subtimestep',
+     &                '','',subtimestep)
+             ! ENDDO !subsurface ice subslope
+
+
 
 c             Starting upward calculations for water :
 c             Actualisation de h2o_vap dans le premier niveau
              zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
-
 c    Take into account the H2O latent heat impact on the surface temperature
               if (latentheat_surfwater) then
@@ -1009 +1131 @@
      &                              *zq_tmp_vap(ig,ilay-1,iq)
               ENDDO
-
 c             Subtimestep water budget :
 
@@ -1016 +1137 @@
               zqsurf(ig)= zqsurf(ig)+(
      &                       zdqsdif(ig))*subtimestep
-
+              if (zqsurf(ig)<0) then 
+                      zqsurf(ig)=0
+              endif
 c             Monitoring instantaneous latent heat flux in W.m-2 :
               zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
@@ -1033 +1156 @@
 c             Fin du sous pas de temps
             ENDDO ! tsub=1,nsubtimestep 
-
 c             Integration of subtimestep temp and water budget :
 c             (btw could also compute the post timestep temp and ice
@@ -1044 +1166 @@
      &                       cos(pi*def_slope_mean(islope)/180.))
      &                       /ptimestep
-
 c if subliming more than qsurf(ice) and on watercaptag, water
 c sublimates from watercap reservoir (dwatercap_dif is <0)
@@ -1063 +1184 @@
            ENDDO ! of DO ig=1,ngrid
           ENDDO ! islope
-
 c Some grid box averages: interface with the atmosphere
        DO ig = 1,ngrid
@@ -1075 +1195 @@
          ENDDO
        ENDDO
-
 ! Recompute values in kg/m^2 slopped
         DO ig = 1,ngrid
@@ -1167 +1286 @@
      &                          "Ground ice latent heat flux",
      &                               "W.m-2",surf_h2o_lh(:,iflat))
-
+                       call write_output('zdq_subtimestep',
+     &                '','',zdqsdif(:)*subtimestep)
+                       call write_output('zdq_end',
+     &                '','',zdqsdif(:))
+                       call write_output('vdifc_subtimestep',
+     &                '','',subtimestep)
 C       Diagnostic output for HDO
 !        if (hdo) then