Ignore:
Timestamp:
Nov 30, 2022, 11:29:29 AM (2 years ago)
Author:
romain.vande
Message:

Mars PEM:
Introduction of the possibility to follow an orbital forcing.
Introduction of new control parameters.
Cleaning of the PEM (removing unused files, add comments and new files)

A file named run_PEM.def can be added to the run.def. It contains the following variables:

_ evol_orbit_pem: Boolean. Do you want to follow an orbital forcing predefined (read in ob_ex_lsp.asc for example)? (default=false)
_ year_bp_ini: Integer. Number of year before present to start the pem run if evol_orbit_pem=.true. , default=0
_ Max_iter_pem: Integer. Maximal number of iteration if none of the stopping criterion is reached and if evol_orbit_pem=.false., default=99999999
_ dt_pem: Integer. Time step of the PEM in year, default=1
_ alpha_criterion: Real. Acceptance rate of sublimating ice surface change, default=0.2
_ soil_pem: Boolean. Do you want to run with subsurface physical processes in the PEM? default=.true.

RV

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2794 r2835  
    11      MODULE ini_soil_mod
    2 
    32
    43      IMPLICIT NONE
    … …  
    65      CONTAINS   
    76
    8 
    97     subroutine ini_icetable(timelen,ngrid,nsoil_PEM, &
    108               therm_i, timestep,tsurf_ave,tsoil_ave,tsurf_inst, tsoil_inst,q_co2,q_h2o,ps,ice_table)
    11 
    12 
    139
    1410      use vertical_layers_mod, only: ap,bp
    … …  
    2824
    2925#include "dimensions.h"
    30 !#include "dimphys.h"
    31 
    32 !#include"comsoil.h"
    33 
    3426
    3527!-----------------------------------------------------------------------
    … …  
    4133      integer,intent(in) :: nsoil_PEM   ! number of soil layers  in the PEM
    4234
    43 
    4435      real,intent(in) :: timestep           ! time step
    4536      real,intent(in) :: tsurf_ave(ngrid)   ! surface temperature
    … …  
    4940      real,intent(in) :: ps(ngrid,timelen)                        ! surface pressure [Pa]
    5041      real,intent(in) :: tsurf_inst(ngrid,timelen) ! soil (mid-layer) temperature
    51 
    5242
    5343! outputs:
    … …  
    5747      real,intent(inout) :: therm_i(ngrid,nsoil_PEM) ! thermal inertia
    5848
    59 
    60      
    6149! local variables:
    6250      integer ig,isoil,it,k,iref
    … …  
    9987    real,allocatable :: diff_rho(:)                    ! difference of vapor content
    10088
    101 
    102       A =(1/m_co2 - 1/m_noco2)
    103       B=1/m_noco2
    104 
    105 
    106  ice_table(:) = 1.
    107 do ig = 1,ngrid
    108  
    109  error_depth = 1.
    110  countloop = 0
    111 
    112 
    113    
    114  do while(( error_depth.gt.tol_error).and.(countloop.lt.countmax).and.(ice_table(ig).gt.-1e-20))
    115    
    116     countloop = countloop +1
    117     Tcol_saved(:) = tsoil_ave(ig,:)
     89    A =(1/m_co2 - 1/m_noco2)
     90    B=1/m_noco2
     91
     92    ice_table(:) = 1.
     93    do ig = 1,ngrid
     94      error_depth = 1.
     95      countloop = 0
     96   
     97      do while(( error_depth.gt.tol_error).and.(countloop.lt.countmax).and.(ice_table(ig).gt.-1e-20))
     98   
     99        countloop = countloop +1
     100        Tcol_saved(:) = tsoil_ave(ig,:)
    118101
    119102!2.  Compute ice table
    120103
    121104! 2.1 Compute  water density at the surface, yearly averaged
    122     allocate(mass_mean(timelen))
     105        allocate(mass_mean(timelen))
    123106!   1.1 Compute the partial pressure of vapor
    124107! a. the molecular mass into the column
    125108       mass_mean(:) = 1/(A*q_co2(ig,:) +B)
    126109! b. pressure level
    127      allocate(zplev(timelen))
    128      do it = 1,timelen
     110        allocate(zplev(timelen))
     111        do it = 1,timelen
    129112         zplev(it) = ap(1) + bp(1)*ps(ig,it)
    130      enddo
     113        enddo
    131114
    132115! c. Vapor pressure
    133      allocate(pvapor(timelen))
    134      pvapor(:) = mass_mean(:)/m_h2o*q_h2o(ig,:)*zplev(:)
    135      deallocate(zplev)
    136      deallocate(mass_mean)
    137 
    138 
     116        allocate(pvapor(timelen))
     117        pvapor(:) = mass_mean(:)/m_h2o*q_h2o(ig,:)*zplev(:)
     118        deallocate(zplev)
     119        deallocate(mass_mean)
    139120 
    140121!  d!  Check if there is frost at the surface and then compute the density
    141122!    at the surface
    142      allocate(rhovapor(timelen))
     123        allocate(rhovapor(timelen))
    143124
    144125         do it = 1,timelen
    … …  
    146127           rhovapor(it) = min(psv_surf,pvapor(it))/tsurf_inst(ig,it)
    147128         enddo
    148      deallocate(pvapor)
    149      rhovapor_avg =     SUM(rhovapor(:),1)/timelen                     
    150      deallocate(rhovapor)
     129        deallocate(pvapor)
     130        rhovapor_avg =SUM(rhovapor(:),1)/timelen                     
     131        deallocate(rhovapor)
    151132
    152133! 2.2 Compute  water density at the soil layer, yearly averaged
    153134
    154      allocate(rho_soil(timelen))
    155      allocate(rho_soil_avg(nsoil_PEM))
    156 
    157 
    158 
    159       do isoil = 1,nsoil_PEM
    160         do it = 1,timelen
    161               rho_soil(it) = exp(alpha/tsoil_inst(ig,isoil,it) +beta)/tsoil_inst(ig,isoil,it)       
    162        enddo
    163        rho_soil_avg(isoil) = SUM(rho_soil(:),1)/timelen     
    164       enddo
    165      deallocate(rho_soil)
     135        allocate(rho_soil(timelen))
     136        allocate(rho_soil_avg(nsoil_PEM))
     137
     138        do isoil = 1,nsoil_PEM
     139          do it = 1,timelen
     140            rho_soil(it) = exp(alpha/tsoil_inst(ig,isoil,it) +beta)/tsoil_inst(ig,isoil,it)       
     141          enddo
     142         rho_soil_avg(isoil) = SUM(rho_soil(:),1)/timelen     
     143        enddo
     144        deallocate(rho_soil)
    166145         
    167 
    168146!2.3 Final: compute ice table
    169          icedepth_prev = ice_table(ig)
    170          ice_table(ig) = -1
    171          allocate(diff_rho(nsoil_PEM))
    172          do isoil = 1,nsoil_PEM
    173              diff_rho(isoil) = rhovapor_avg -  rho_soil_avg(isoil)
     147        icedepth_prev = ice_table(ig)
     148        ice_table(ig) = -1
     149        allocate(diff_rho(nsoil_PEM))
     150        do isoil = 1,nsoil_PEM
     151          diff_rho(isoil) = rhovapor_avg -  rho_soil_avg(isoil)
     152        enddo
     153        deallocate(rho_soil_avg)
     154
     155        if(diff_rho(1) > 0) then
     156          ice_table(ig) = 0.
     157        else
     158          do isoil = 1,nsoil_PEM -1
     159            if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
     160              z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
     161              z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
     162              ice_table(ig) = -z2/z1
     163              exit
     164            endif
     165          enddo
     166        endif
     167   
     168        deallocate(diff_rho)
     169
     170!3. Update Soil Thermal Inertia
     171
     172        if (ice_table(ig).gt. 0.) then
     173          if (ice_table(ig).lt. 1e-10) then
     174            do isoil = 1,nsoil_PEM
     175              therm_i(ig,isoil)=ice_inertia
     176            enddo
     177          else
     178      ! 4.1 find the index of the mixed layer
     179          iref=0 ! initialize iref
     180          do k=1,nsoil_PEM ! loop on layers
     181            if (ice_table(ig).ge.layer_PEM(k)) then
     182              iref=k ! pure regolith layer up to here
     183            else
     184              ! correct iref was obtained in previous cycle
     185            exit
     186            endif
     187          enddo
     188
     189      ! 4.2 Build the new ti
     190         do isoil=1,iref
     191           therm_i(ig,isoil) =inertiedat_PEM(ig,isoil)
    174192         enddo
    175          deallocate(rho_soil_avg)
    176 
    177 
    178          if(diff_rho(1) > 0) then
    179            ice_table(ig) = 0.
    180          else
    181            do isoil = 1,nsoil_PEM -1
    182              if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
    183                z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
    184                z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
    185                ice_table(ig) = -z2/z1
    186                exit
    187              endif
    188            enddo
    189           endif
    190 
    191 
    192    
    193     deallocate(diff_rho)
    194 
    195 
    196 !3. Update Soil Thermal Inertia
    197 
    198   if (ice_table(ig).gt. 0.) then
    199   if (ice_table(ig).lt. 1e-10) then
    200       do isoil = 1,nsoil_PEM
    201       therm_i(ig,isoil)=ice_inertia
    202       enddo
    203   else
    204       ! 4.1 find the index of the mixed layer
    205       iref=0 ! initialize iref
    206       do k=1,nsoil_PEM ! loop on layers
    207         if (ice_table(ig).ge.layer_PEM(k)) then
    208           iref=k ! pure regolith layer up to here
    209         else
    210          ! correct iref was obtained in previous cycle
    211          exit
    212         endif
    213        
    214       enddo
    215 
    216 
    217 
    218       ! 4.2 Build the new ti
    219       do isoil=1,iref
    220          therm_i(ig,isoil) =inertiedat_PEM(ig,isoil)
    221       enddo
    222 
    223 
    224       if (iref.lt.nsoil_PEM) then
    225         if (iref.ne.0) then
    226           ! mixed layer
    227            therm_i(ig,iref+1)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
    228             (((ice_table(ig)-layer_PEM(iref))/(inertiedat_PEM(ig,iref)**2))+ &
     193
     194         if (iref.lt.nsoil_PEM) then
     195           if (iref.ne.0) then
     196             ! mixed layer
     197             therm_i(ig,iref+1)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
     198              (((ice_table(ig)-layer_PEM(iref))/(inertiedat_PEM(ig,iref)**2))+ &
    229199                      ((layer_PEM(iref+1)-ice_table(ig))/(ice_inertia**2))))
    230200           
    231 
    232 
    233         else ! first layer is already a mixed layer
    234           ! (ie: take layer(iref=0)=0)
    235           therm_i(ig,1)=sqrt((layer_PEM(1))/ &
     201           else ! first layer is already a mixed layer
     202            ! (ie: take layer(iref=0)=0)
     203            therm_i(ig,1)=sqrt((layer_PEM(1))/ &
    236204                          (((ice_table(ig))/(inertiedat_PEM(ig,1)**2))+ &
    237205                           ((layer_PEM(1)-ice_table(ig))/(ice_inertia**2))))
    238         endif ! of if (iref.ne.0)       
    239         ! lower layers of pure ice
    240         do isoil=iref+2,nsoil_PEM
    241           therm_i(ig,isoil)=ice_inertia
    242         enddo
    243       endif ! of if (iref.lt.(nsoilmx))
    244       endif ! permanent glaciers
    245 
    246 
     206           endif ! of if (iref.ne.0)       
     207           ! lower layers of pure ice
     208           do isoil=iref+2,nsoil_PEM
     209             therm_i(ig,isoil)=ice_inertia
     210           enddo
     211         endif ! of if (iref.lt.(nsoilmx))
     212       endif ! permanent glaciers
    247213
    248214       call soil_pem_1D(nsoil_PEM,.true.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM)
    249 
    250215       call soil_pem_1D(nsoil_PEM,.false.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM)
    251216
    252 
    253       do it = 1,timelen
    254         tsoil_inst(ig,:,it) = tsoil_inst(ig,:,it) - (Tcol_saved(:) - tsoil_ave(ig,:))
    255       enddo
    256 
    257 
    258 
    259       error_depth = abs(icedepth_prev - ice_table(ig))
    260      
    261      endif
    262 
    263 
    264 enddo
    265 
    266  error_depth = 1.
    267  countloop = 0
    268 enddo
     217       do it = 1,timelen
     218         tsoil_inst(ig,:,it) = tsoil_inst(ig,:,it) - (Tcol_saved(:) - tsoil_ave(ig,:))
     219       enddo
     220
     221       error_depth = abs(icedepth_prev - ice_table(ig))     
     222       endif
     223
     224     enddo
     225
     226     error_depth = 1.
     227     countloop = 0
     228   enddo
    269229 
    270 
    271230      END SUBROUTINE ini_icetable
    272231      subroutine soil_pem_1D(nsoil,firstcall, &
    273232               therm_i,                          &
    274233               timestep,tsurf,tsoil,alph,beta)
    275 
    276234
    277235      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
    … …  
    309267      real,intent(inout) :: alph(nsoil-1)
    310268      real,intent(inout) :: beta(nsoil-1)
    311 
    312 
    313 
    314 
    315269     
    316270! local variables:
    … …  
    330284        enddo
    331285
    332 
    333286! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
    334287        do ik=1,nsoil-1
    … …  
    336289                     +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ik-1))  &
    337290                         /(mlayer_PEM(ik)-mlayer_PEM(ik-1))
    338 
    339291        enddo
    340292
    … …  
    365317        enddo
    366318
    367 
    368 
    369      
    370  
    371319endif ! of if (firstcall)
    372 
    373 
    374320
    375321      IF (.not.firstcall) THEN
    … …  
    387333     
    388334          ENDIF
    389 
    390 
    391 
    392335
    393336!  2. Compute beta_PEM coefficients (preprocessing for next time step)
    … …  
    405348      enddo
    406349
    407 
    408 
    409350      end
    410351
    411 
    412 
    413 
    414       END MODULE ini_soil_mod
    415 
    416 
    417      
     352      END MODULE ini_soil_mod     
Note: See TracChangeset for help on using the changeset viewer.