Changeset 3123 for trunk/LMDZ.COMMON


Ignore:
Timestamp:
Nov 11, 2023, 5:34:53 PM (13 months ago)
Author:
llange
Message:

MARS PEM
1) Adapting the Mars PEM soil grid to the one of the PCM. The first layers in the PEM follow those from the PCM (PYM grid), and then, for layers at depth, we use the classical power low grid.
2) Correction when reading the soil temperature, water density at the surface, and initialization of the ice table.
LL

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3122 r3123  
    133133== 10/11/2023 == JBC
    134134Correction of the reading of the PCM data (it did not work if no slope was used) + some minor related cleanings.
     135
     136== 11/11/2023 == LL
     137Adapting the PEM soil grid to the one of the PCM
     138Minor corrections when reading/initializing soil temperature, subsurface water ice

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

    r2980 r3123  
    22
    33implicit none
    4   integer, parameter :: nsoilmx_PEM = 27                 ! number of layers in the PEM
     4  integer, parameter :: nsoilmx_PEM = 68                 ! number of layers in the PEM
    55  real,save,allocatable,dimension(:) :: layer_PEM        ! soil layer depths   [m]
    66  real,save,allocatable,dimension(:) :: mlayer_PEM       ! soil mid-layer depths [m]

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

    r3122 r3123  
    11371137             co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
    11381138write(*,*) "restartpem.nc has been written"
    1139 
    11401139call info_PEM(year_iter,criterion_stop,i_myear,n_myear)
    11411140

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

    r3082 r3123  
    262262      write(*,*)'PEM settings: failed loading <ice_table>'
    263263      write(*,*)'will reconstruct the values of the ice table given the current state'
     264      ice_table(:,:) = -1  ! by default, no ice table 
     265      ice_table_thickness(:,:) = -1  ! by default, no ice table
    264266     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, TI_PEM(:,1,:),ice_table,ice_table_thickness)
    265267     call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
    … …  
    427429      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    428430      call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    429 
    430 
    431     do it = 1,timelen
     431 
     432! First raw initialization
     433      do it = 1,timelen
    432434        do isoil = nsoil_GCM+1,nsoil_PEM
    433       call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_inst(:,:,islope,it))
     435          tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
    434436        enddo
    435      enddo
    436 
     437      enddo
     438
     439      do it = 1,timelen
    437440        do isoil = nsoil_GCM+1,nsoil_PEM
    438           do ig = 1,ngrid
    439             watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
    440           enddo
     441          call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_inst(:,:,islope,it))
    441442        enddo
     443       enddo
     444
     445       do isoil = nsoil_GCM+1,nsoil_PEM
     446         do ig = 1,ngrid
     447           watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
     448         enddo
     449       enddo
    442450enddo !islope
    443451write(*,*) 'PEMETAT0: TSOIL done'
    … …  
    445453!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    446454!c) Ice table
     455       ice_table(:,:) = -1.  ! by default, no ice table
     456       ice_table_thickness(:,:) = -1.
    447457       call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
    448458       call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)

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

    r3122 r3123  
    6464integer                         :: islope              ! loop for variables
    6565character(2)                    :: num                 ! for reading sloped variables
    66 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn         ! h2o ice per slope of the concatenated file [kg/m^2]
     66real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn             ! h2o ice per slope of the concatenated file [kg/m^2]
    6767real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
    68 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm           ! CO2 volume mixing ratio in the first layer [mol/m^3]
    69 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                ! Surface Pressure [Pa]
     68real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm               ! CO2 volume mixing ratio in the first layer [mol/m^3]
     69real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                    ! Surface Pressure [Pa]
    7070real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
    7171real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
    72 real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn         ! Average surface temperature of the concatenated file [K]
    73 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn         ! Average soil temperature of the concatenated file [K]
    74 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn         ! Surface temperature of the concatenated file, time series [K]
    75 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn         ! Soil temperature of the concatenated file, time series [K]
    76 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn             ! CO2 mass mixing ratio in the first layer [kg/m^3]
    77 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn             ! H2O mass mixing ratio in the first layer [kg/m^3]
    78 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn     ! co2 ice amount per  slope of the year [kg/m^2]
    79 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn ! Water density at the surface, time series [kg/m^3]
    80 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn ! Water density in the soil layer, time series [kg/m^3]
    81 real, dimension(ngrid,nslope,timelen)                               :: watersurf_density     ! Water density at the surface, time series [kg/m^3]
     72real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn             ! Average surface temperature of the concatenated file [K]
     73real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn             ! Average soil temperature of the concatenated file [K]
     74real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn             ! Surface temperature of the concatenated file, time series [K]
     75real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn             ! Soil temperature of the concatenated file, time series [K]
     76real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn                 ! CO2 mass mixing ratio in the first layer [kg/m^3]
     77real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn                 ! H2O mass mixing ratio in the first layer [kg/m^3]
     78real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn         ! co2 ice amount per  slope of the year [kg/m^2]
     79real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn     ! Water density at the surface, time series [kg/m^3]
     80real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: watersurf_density_dyn_ave ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3]
     81real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn     ! Water density in the soil layer, time series [kg/m^3]
     82real, dimension(ngrid,nslope,timelen)                               :: watersurf_density         ! Water density at the surface, physical grid, time series [kg/m^3]
    8283!-----------------------------------------------------------------------
    8384! Open the NetCDF file
    … …  
    115116        write(*,*) "Data for soiltemp downloaded!"
    116117
    117         call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
     118        call get_var4("Waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
    118119        write(*,*) "Data for waterdensity_soil downloaded!"
    119120
    120         call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,islope,:))
     121        call get_var3("Waterdensity_surface",watersurf_density_dyn(:,:,1,:))
    121122        write(*,*) "Data for waterdensity_surface downloaded!"
    122123    endif !soil_pem
    … …  
    184185tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
    185186
    186 #ifndef CPP_1D
    187     do islope = 1,nslope
    188         do t = 1,timelen
    189             call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
    190         enddo
    191     enddo
    192 #endif
    193 
    194187if (soil_pem) then
    195188    write(*,*) "Computing average of tsoil..."
    196189    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
    197190    write(*,*) "Computing average of waterdensity_surface..."
    198     watersurf_density_ave(:,:) = sum(watersurf_density(:,:,:),3)/timelen
     191    watersurf_density_dyn_ave(:,:,:) = sum(watersurf_density_dyn(:,:,:,:),4)/timelen
    199192endif
    200193
    … …  
    238231                enddo
    239232            enddo
     233            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_ave(:,:,islope),watersurf_density_ave(:,islope))
    240234        endif ! soil_pem
    241235        do t = 1,timelen
    … …  
    257251        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
    258252        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
     253        watersurf_density_ave(1,:) = watersurf_density_dyn_ave(1,1,:)
    259254    endif ! soil_pem
    260255    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)

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

    r3076 r3123  
    4949real    :: alpha,lay1 ! coefficients for building layers
    5050real    :: xmin, xmax ! to display min and max of a field
     51real    :: index_powerlaw ! coefficient to match the power low grid with the exponential one
    5152!=======================================================================
    5253! 1. Depth coordinate
    5354! -------------------
    54 ! 1.4 Build mlayer(), if necessary
    55 ! default mlayer distribution, following a power law:
    56 !  mlayer(k)=lay1*alpha**(k-1/2)
     55! mlayer_PEM distribution: For the grid in common with the PCM: lay1*(1+k^(2.9)*(1-exp(-k/20))
     56! Then we use a power low : lay1*alpha**(k-1/2)
    5757lay1 = 2.e-4
    5858alpha = 2
    59 do iloop = 0,nsoil_PEM - 1
    60     mlayer_PEM(iloop) = lay1*(alpha**(iloop - 0.5))
     59do iloop = 0,nsoil_GCM - 1
     60    mlayer_PEM(iloop) = lay1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
    6161enddo
    6262
    63 ! 1.5 Build layer(); following the same law as mlayer()
    64 ! Assuming layer distribution follows mid-layer law:
    65 ! layer(k)=lay1*alpha**(k-1)
    66 lay1 = sqrt(mlayer_PEM(0)*mlayer_PEM(1))
    67 alpha = mlayer_PEM(1)/mlayer_PEM(0)
    68 do iloop = 1,nsoil_PEM
    69     layer_PEM(iloop) = lay1*(alpha**(iloop - 1))
     63do iloop = 0, nsoil_PEM-1
     64    if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_GCM-1)) then
     65        index_powerlaw = iloop
     66        exit
     67    endif
    7068enddo
     69
     70do iloop = nsoil_GCM, nsoil_PEM-1
     71    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop-nsoil_GCM)-0.5))
     72enddo
     73
     74! Build layer_PEM()
     75do iloop=1,nsoil_PEM-1
     76layer_PEM(iloop)=(mlayer_PEM(iloop)+mlayer_PEM(iloop-1))/2.
     77enddo
     78layer_PEM(nsoil_PEM)=2*mlayer_PEM(nsoil_PEM-1) - mlayer_PEM(nsoil_PEM-2)
     79
    7180
    7281! 2. Thermal inertia (note: it is declared in comsoil_h)
Note: See TracChangeset for help on using the changeset viewer.