Changeset 2902


Ignore:
Timestamp:
Feb 22, 2023, 6:56:01 PM (21 months ago)
Author:
llange
Message:

PEM
First implementation of a dynamic ice table following Schorgofer, Icarus (2010) / MSIM
The dynamic of the ice is implemented here. Impact on the soil properties (thermal, inertia) is currently implemented and will be commit later
Small fix for the ice table at equilibrium
LL

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90

    r2888 r2902  
    2828
    2929#ifndef CPP_STD
    30     USE comsoil_h_PEM, only: layer_PEM                             ! Depth of the vertical grid
     30    USE comsoil_h_PEM, only: mlayer_PEM                             ! Depth of the vertical grid
    3131    implicit none
    3232
    … …  
    4646      else
    4747        do islope = 1,nslope
    48            ice_table(ig,islope) = -10.
     48           ice_table(ig,islope) = -1.
    4949
    5050         do isoil = 1,nsoil_PEM
    … …  
    5757           do isoil = 1,nsoil_PEM -1                                ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
    5858             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
    59                z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
    60                z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
    61                ice_table(ig,islope) = -z2/z1
     59               call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table(ig,islope))
    6260               exit
    6361             endif !diff_rho(z) <0 & diff_rho(z+1) > 0
    … …  
    7270      END
    7371
    74 
    75 
    76 end module
     72   SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
     73!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     74!!!
     75!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
     76!!!
     77!!! Author: LL
     78!!!
     79!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     80use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM
     81implicit none
     82! inputs
     83! ------
     84
     85real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
     86real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa]
     87real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa]
     88real,intent(in) :: pwat_surf ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
     89real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
     90real,intent(in) :: B ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
     91integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1]
     92real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
     93
     94! outputs
     95! -------
     96integer,intent(out) :: index_filling ! index where the pore filling begins [1]
     97integer, intent(out) :: index_geothermal ! index where the ice table stops [1]
     98real, intent(out) :: depth_geothermal ! depth where the ice table stops [m]
     99real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
     100
     101! local
     102
     103real :: eta(nsoilmx_PEM)  ! constriction
     104integer :: ilay ! index for loop
     105real :: old_depth_filling ! depth_filling saved
     106real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn
     107integer :: index_tmp ! for loop
     108real :: Jdry ! flux trought the dry layer
     109real :: Jsat ! flux trought the ice layer
     110real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer
     111integer :: index_firstice ! first index where ice appears (i.e., f > 0)
     112real :: dz_eta(nsoilmx_PEM) ! \partial z \eta
     113real :: dz_eta_low ! same but evaluated at the interface for ice
     114real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat
     115real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal
     116
     117! constant
     118real :: pvap2rho = 18.e-3/8.314
     119!
     120
     121
     122
     123
     124 
     125! 0. Compute constriction over the layer
     126! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
     127if (index_IS.lt.0) then
     128   index_tmp = nsoilmx_PEM
     129   do ilay = 1,nsoilmx_PEM
     130       call constriction(porefill(ilay),eta(ilay))
     131   enddo
     132else
     133   index_tmp = index_IS
     134   do ilay = 1,index_IS-1
     135       call constriction(porefill(ilay),eta(ilay))
     136   enddo
     137   do ilay = index_IS,nsoilmx_PEM
     138       eta(ilay) = 0.
     139   enddo
     140endif
     141
     142! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)   
     143
     144old_depth_filling = depth_filling
     145
     146call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) 
     147
     148do ilay = 1,index_tmp
     149  Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
     150  Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
     151  if((Jdry - Jsat).le.0) then
     152     index_filling = ilay
     153     exit
     154   endif
     155enddo
     156
     157if(index_filling.eq.1)  depth_filling = mlayer_PEM(1)
     158
     159if(index_filling.gt.1) then
     160    Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1)
     161    Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1)
     162   call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling)
     163endif
     164
     165
     166! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
     167
     168
     169! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
     170
     171index_firstice = -1
     172do ilay = 1,nsoilmx_PEM
     173   if (porefill(ilay).le.0.) then
     174        index_firstice = ilay  ! first point with ice
     175        exit
     176    endif
     177enddo
     178
     179! 2.1: now we can computeCompute d_z (eta* d_z(rho))
     180
     181call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta)
     182
     183if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then
     184     call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
     185endif
     186
     187call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
     188dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:))
     189
     190! 3. Ice table boundary due to geothermal heating
     191
     192if(index_IS.gt.0) index_geothermal = -1
     193if(index_geothermal.lt.0) depth_geothermal = -1.
     194if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
     195   index_geothermal = -1
     196   do ilay=2,nsoilmx_PEM
     197        if (dz_psat(ilay).gt.0.) then  ! first point with reversed flux
     198           index_geothermal=ilay
     199           call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
     200           exit
     201        endif
     202     enddo
     203  else
     204     index_geothermal = -1
     205  endif
     206 if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then !  Eq. 24 from Schorgofer, Icarus (2010)
     207     call  colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove)
     208     index_tmp = -1
     209     do ilay=index_geothermal,nsoilmx_PEM
     210        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle  ! eta=0 means completely full
     211         call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
     212        if (massfillafter<dz_psat(ilay)*pvap2rho*B) then  ! usually executes on i=typeG
     213           if (ilay>index_geothermal) then
     214!              write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
     215              call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, & 
     216                   dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
     217!               if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*)
     218!                     '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay)
     219              index_tmp=ilay
     220           endif
     221           ! otherwise leave depth_geothermal unchanged
     222           exit
     223        endif
     224        massfillabove = massfillafter
     225     enddo
     226     if (index_tmp.gt.0) index_geothermal = index_tmp
     227  end if
     228  return
     229end subroutine
     230
     231
     232
     233
     234
     235SUBROUTINE findroot(y1,y2,z1,z2,zr)
     236        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     237        !!!
     238        !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
     239        !!!
     240        !!! Author: LL
     241        !!!
     242        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     243
     244        implicit none
     245        real,intent(in) :: y1,y2 ! difference between surface water density and at depth  [kg/m^3]
     246        real,intent(in) :: z1,z2 ! depth [m]
     247        real,intent(out) :: zr ! depth at which we have zero
     248        zr = (y1*z2 - y2*z1)/(y1-y2)
     249        RETURN
     250        end
     251
     252SUBROUTINE constriction(porefill,eta)
     253
     254        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     255        !!!
     256        !!! Purpose: Compute the  constriction of vapor flux by pore ice
     257        !!!
     258        !!! Author: LL
     259        !!!
     260        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     261        implicit none
     262        real,intent(in) :: porefill ! pore filling fraction
     263        real,intent(out) :: eta ! constriction
     264
     265        !!!
     266
     267
     268        if (porefill.le.0.) eta = 1.
     269        if ((porefill.gt.0.) .and.(porefill.lt.1.)) then
     270        eta = (1-porefill)**2  ! Hudson et al., JGR, 2009
     271        endif
     272        if (porefill.le.1.) eta = 0.
     273        return
     274        end
     275
     276
     277
     278subroutine deriv1(z,nz,y,y0,ybot,dzY)
     279!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     280!!!
     281!!! Purpose: Compute the first derivative of a function y(z) on an irregular grid
     282!!!           
     283!!! Author: From N.S (N.S, Icarus 2010), impletented here by LL
     284!!!
     285!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     286        implicit none
     287
     288! first derivative of a function y(z) on irregular grid
     289! upper boundary conditions: y(0)=y0
     290! lower boundary condition.: yp = ybottom
     291  integer, intent(IN) :: nz ! number of layer
     292  real, intent(IN) :: z(nz) ! depth layer
     293  real, intent(IN) :: y(nz) ! function which needs to be derived
     294  real, intent(IN) :: y0,ybot ! boundary conditions
     295  real, intent(OUT) :: dzY(nz) ! derivative of y w.r.t depth
     296! local
     297  integer :: j
     298  real :: hm,hp,c1,c2,c3
     299
     300  hp = z(2)-z(1)
     301  c1 = z(1)/(hp*z(2))
     302  c2 = 1/z(1) - 1/(z(2)-z(1))
     303  c3 = -hp/(z(1)*z(2))
     304  dzY(1) = c1*y(2) + c2*y(1) + c3*y0
     305  do j=2,nz-1
     306     hp = z(j+1)-z(j)
     307     hm = z(j)-z(j-1)
     308     c1 = +hm/(hp*(z(j+1)-z(j-1)))
     309     c2 = 1/hm - 1/hp
     310     c3 = -hp/(hm*(z(j+1)-z(j-1)))
     311     dzY(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
     312  enddo
     313  dzY(nz) = (ybot - y(nz-1))/(2.*(z(nz)-z(nz-1)))
     314 return
     315end subroutine deriv1
     316
     317
     318
     319subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
     320!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     321!!!
     322!!! Purpose: Compute the second derivative of a function y(z) on an irregular grid
     323!!!           
     324!!! Author: N.S (raw copy/paste from MSIM)
     325!!!
     326!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     327
     328  ! second derivative y_zz on irregular grid
     329  ! boundary conditions: y(0)=y0, y(nz+1)=yNp1
     330  implicit none
     331  integer, intent(IN) :: nz
     332  real, intent(IN) :: z(nz),y(nz),y0,yNp1
     333  real, intent(OUT) :: yp2(nz)
     334  integer j
     335  real hm,hp,c1,c2,c3
     336
     337  c1 = +2./((z(2)-z(1))*z(2))
     338  c2 = -2./((z(2)-z(1))*z(1))
     339  c3 = +2./(z(1)*z(2))
     340  yp2(1) = c1*y(2) + c2*y(1) + c3*y0
     341
     342  do j=2,nz-1
     343     hp = z(j+1)-z(j)
     344     hm = z(j)-z(j-1)
     345     c1 = +2./(hp*(z(j+1)-z(j-1)))
     346     c2 = -2./(hp*hm)
     347     c3 = +2./(hm*(z(j+1)-z(j-1)))
     348     yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
     349  enddo
     350  yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2
     351  return
     352end subroutine deriv2_simple
     353
     354
     355
     356subroutine  deriv1_onesided(j,z,nz,y,dy_zj)
     357!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     358!!!
     359!!! Purpose:   First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
     360!!!           
     361!!! Author: N.S (raw copy/paste from MSIM)
     362!!!
     363!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     364
     365  implicit none
     366  integer, intent(IN) :: nz,j
     367  real, intent(IN) :: z(nz),y(nz)
     368  real, intent(out) :: dy_zj
     369  real h1,h2,c1,c2,c3
     370  if (j<1 .or. j>nz-2) then
     371     dy_zj = -1.
     372  else
     373     h1 = z(j+1)-z(j)
     374     h2 = z(j+2)-z(j+1)
     375     c1 = -(2*h1+h2)/(h1*(h1+h2))
     376     c2 =  (h1+h2)/(h1*h2)
     377     c3 = -h1/(h2*(h1+h2))
     378     dy_zj = c1*y(j) + c2*y(j+1) + c3*y(j+2)
     379  endif
     380  return
     381end subroutine deriv1_onesided
     382
     383
     384subroutine  colint(y,z,nz,i1,i2,integral)
     385!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     386!!!
     387!!! Purpose:  Column integrates y on irregular grid
     388!!!           
     389!!! Author: N.S (raw copy/paste from MSIM)
     390!!!
     391!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     392 
     393  implicit none
     394  integer, intent(IN) :: nz, i1, i2
     395  real, intent(IN) :: y(nz), z(nz)
     396  real,intent(out) :: integral
     397  integer i
     398  real dz(nz)
     399 
     400  dz(1) = (z(2)-0.)/2
     401  do i=2,nz-1
     402     dz(i) = (z(i+1)-z(i-1))/2.
     403  enddo
     404  dz(nz) = z(nz)-z(nz-1)
     405  integral = sum(y(i1:i2)*dz(i1:i2))
     406end subroutine colint
     407
     408
     409
     410
     411
     412
     413        end module

  • trunk/LMDZ.MARS/README

    r2900 r2902  
    39043904nslope is now read and written in the startfi.nc
    39053905
    3906 
     3906== 22/02/2023 == LL
     3907PEM: First implementation of a dynamic ice table (vs. static ice table at equilibrium) following Schorgofer (2010,Icarus)

  • trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F

    r2823 r2902  
    170170      REAL masseq(nlayer),wq(nlayer+1)
    171171      INTEGER ifils,iq2
    172 
     172c LL: flux dependant albedo
     173      real :: totflux(ngrid)
    173174c----------------------------------------------------------------------
    174175
    … …  
    176177c   --------------
    177178c
     179       totflux(:) = fluxsurf_sw(:,1)+fluxsurf_sw(:,2)
    178180      ! AS: firstcall OK absolute
    179181      IF (firstcall) THEN
    … …  
    500502!     ----------------------------------------
    501503      CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
     504      DO ig = 1,ngrid
     505        if(piceco2(ig).gt.0.) then
     506          psolaralb(ig,1) = 0.267+0.00059*totflux(ig)
     507          psolaralb(ig,2) = psolaralb(ig,1)
     508          if(ig.eq.ngrid) then
     509           psolaralb(ig,1) = 0.532+8.72e-4*totflux(ig)
     510           psolaralb(ig,2) = psolaralb(ig,1)
     511          endif
     512        endif
     513      ENDDO
     514
     515
    502516
    503517! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
    … …  
    739753!             write(*,*) "co2condens: South pole: latitude(ngrid)=",
    740754!     &                                           latitude(ngrid)
    741              ztcondsol(ngrid)=
    742      &          1./(bcond-acond*log(.01*vmr_co2(ngrid,1)*
    743      &                    (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
    744              pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
     755!             ztcondsol(ngrid)=
     756!     &          1./(bcond-acond*log(.01*vmr_co2(ngrid,1)*
     757!     &                    (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
     758!             pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
    745759           endif
    746760         endif
Note: See TracChangeset for help on using the changeset viewer.