Ignore:
Timestamp:
Dec 20, 2022, 7:07:27 PM (2 years ago)
Author:
llange
Message:

PEM
Update the codes for subsurface evolution to run with XIOS
1) Water density at the surface and in the soil is now read in the XIOS file
2) Reshape routine created as XIOS output have one element less in the longitude grid
3) The ice table is now computed using these water densities
4) Minor fixs in the main, adsorption module, and tendencies evolutions

LL

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

Legend:

Unmodified
Added
Removed

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

    r2842 r2849  
    3232 REAL ::  rho_regolith = 2000.    ! density of the reoglith, Buhler & Piqueux 2021
    3333 real :: alpha_clapeyron = -6143.7! eq. (2) in Murphy & Koop 2005
    34  real :: beta_clapeyron = 29.9074 ! eq. (2) in Murphy & Koop 2005
     34 real :: beta_clapeyron = 28.9074 ! eq. (2) in Murphy & Koop 2005
    3535 real :: mi = 2.84e-7             ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
    3636 real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
    … …  
    8888  do islope = 1,nslope
    8989    do iloop = 1,n_1km
     90       K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
    9091     if(TI_PEM(ig,iloop,islope).lt.inertie_thresold)  then
    91        K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
     92
    9293       theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu
    9394       m_h2o_adsorbed(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith

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

    r2842 r2849  
    1    SUBROUTINE computeice_table(timelen,ngrid,nslope,nsoil_GCM,nsoil_PEM,tsoil,tsurf,q_co2,q_h2o,ps,ice_table)
     1   SUBROUTINE computeice_table(ngrid,nslope,nsoil_PEM,rhowatersurf_ave,rhowatersoil_ave,ice_table)
    22#ifndef CPP_STD
    33    USE comsoil_h, only:  inertiedat, volcapa
    … …  
    66    implicit none
    77
    8     integer,intent(in) :: timelen,ngrid,nslope,nsoil_PEM,nsoil_GCM
    9     real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope,timelen)    ! soil temperature, timeseries [K]
    10     real,intent(in) :: tsurf(ngrid,nslope,timelen)              ! surface temperature [K]
    11     real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
    12     real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
    13     real,intent(in) :: ps(ngrid,timelen)                        ! surface pressure [Pa]
    14     real,intent(out) :: ice_table(ngrid,nslope)                 ! ice table [m]
     8    integer,intent(in) :: ngrid,nslope,nsoil_PEM
     9    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)    ! surface temperature, timeseries [K]
     10    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! soil water density, yearly average [kg/m^3]
     11    real,intent(inout) :: ice_table(ngrid,nslope)                 ! ice table [m]
     12    real :: z1,z2
     13    integer ig, islope,isoil
     14    real :: diff_rho(nsoil_PEM)                    ! difference of vapor content
    1515
    16     real :: m_h2o = 18.01528E-3
    17     real :: m_co2 = 44.01E-3 
    18     real :: m_noco2 = 33.37E-3 
    19     real :: A,B,z1,z2
    20     real :: alpha = -6143.7
    21     real :: beta = 29.9074
    22 
    23     integer ig, islope,isoil,it
    24     real,allocatable :: mass_mean(:,:)                            ! mean mass above the surface
    25     real,allocatable :: zplev_mean(:,:)                           ! pressure above the surface
    26     real,allocatable :: pvapor(:,:)                               ! partial pressure above the surface
    27 
    28     real,allocatable :: rhovapor(:,:,:)
    29     real,allocatable :: rhovapor_avg(:,:)                          ! mean vapor_density at the surface yearly averaged
    30 
    31     real :: psv_surf
    32     real,allocatable :: rho_soil(:,:,:,:)            ! water vapor in the soil
    33     real,allocatable :: rho_soil_avg(:,:,:)                ! water vapor in the soil yearly averaged
    34 
    35     real,allocatable :: diff_rho(:,:,:)                    ! difference of vapor content
    36 
    37      allocate(rhovapor(ngrid,nslope,timelen))
    38      allocate(rhovapor_avg(ngrid,nslope))
    39      allocate(pvapor(ngrid,timelen))
    40 
    41      allocate(mass_mean(ngrid,timelen))
    42      allocate(zplev_mean(ngrid,timelen))
    43 
    44 ! 0. Some initializations
    45 
    46       A =(1/m_co2 - 1/m_noco2)
    47       B=1/m_noco2
    48 ! 1. Compute rho surface yearly averaged
    49 
    50 !   1.1 Compute the partial pressure of vapor
    51 !a. the molecular mass into the column
    52      do ig = 1,ngrid
    53        mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
    54      enddo
    55 
    56 ! b. pressure level
    57      do it = 1,timelen
    58        do ig = 1,ngrid
    59          zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
    60        enddo
    61      enddo
    62 
    63 ! c. Vapor pressure
    64      pvapor(:,:) = mass_mean(:,:)/m_h2o*q_h2o(:,:)*zplev_mean(:,:)
    65    
    66      deallocate(mass_mean)
    67      deallocate(zplev_mean)
    68  
    69 !    1.2   Check if there is frost at the surface and then compute the density
    70 !    at the surface
    71      do ig = 1,ngrid
    72        do islope = 1,nslope
    73          do it = 1,timelen
    74            psv_surf = exp(alpha/tsurf(ig,islope,it) +beta)
    75           if ((isnan(psv_surf)).or.(isnan(pvapor(ig,it)))) then
    76           write(*,*) 'pb vapor',ig,islope,it
    77          stop
    78          endif       
    79          rhovapor(ig,islope,it) = min(psv_surf,pvapor(ig,it))/tsurf(ig,islope,it)
    80          enddo
    81        enddo
    82      enddo
    83      deallocate(pvapor)
    84 
    85 !    1.3 Density at the surface yearly averaged
    86      rhovapor_avg(:,:) = SUM(rhovapor(:,:,:),3)/timelen
    87 
    88      deallocate(rhovapor)
    89 
    90 ! 2. Compute rho soil vapor
    91    
    92      allocate(rho_soil_avg(ngrid,nslope,nsoil_PEM))
    93      allocate(rho_soil(ngrid,nslope,nsoil_PEM,timelen))
    9416
    9517     do ig = 1,ngrid
    96        do islope = 1,nslope
    97          do isoil = 1,nsoil_PEM
    98             do it = 1,timelen
    99               rho_soil(ig,islope,isoil,it) = exp(alpha/tsoil(ig,isoil,islope,it) +beta)/tsoil(ig,isoil,islope,it)       
    100              enddo
    101            enddo
    102        enddo
    103      enddo
    104 
    105     rho_soil_avg(:,:,:) = SUM( rho_soil(:,:,:,:),4)/timelen
    106     deallocate(rho_soil)
    107 
    108 ! 3. Computing ice table
    109    
    110     ice_table (:,:) = -1e4
    111 
    112          allocate(diff_rho(ngrid,nslope,nsoil_PEM))
     18      do islope = 1,nslope
     19           ice_table(ig,islope) = -10.
    11320
    11421         do isoil = 1,nsoil_PEM
    115            diff_rho(:,:,isoil) = rhovapor_avg(:,:) - rho_soil_avg(:,:,isoil)
    116 !          write(*,*) 'diff =',ig,islope,isoil,diff_rho(ig,islope,isoil),rhovapor_avg(ig,islope) ,rho_soil_avg(ig,islope,isoil)
     22           diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
     23
    11724         enddo
    118    
    119   deallocate(rhovapor_avg)
    120   deallocate(rho_soil_avg)
    12125
    122      do ig = 1,ngrid
    123        do islope = 1,nslope
    124          if(diff_rho(ig,islope,1) > 0) then
     26         if(diff_rho(1) > 0) then
    12527           ice_table(ig,islope) = 0.
    12628         else
    12729           do isoil = 1,nsoil_PEM -1
    128              if((diff_rho(ig,islope,isoil).lt.0).and.(diff_rho(ig,islope,isoil+1).gt.0.)) then
    129                z1 = (diff_rho(ig,islope,isoil) - diff_rho(ig,islope,isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
    130                z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(ig,islope,isoil+1)
     30             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
     31               z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
     32               z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
    13133               ice_table(ig,islope) = -z2/z1
    13234               exit
    … …  
    13638        enddo
    13739      enddo
     40
    13841   
    139     deallocate(diff_rho)
    14042
    14143!=======================================================================

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

    r2835 r2849  
    3737
    3838!=======================================================================
     39
     40  STOPPING=.false.
    3941
    4042  pos_tend=0.

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

    r2842 r2849  
    7878                           co2iceflow, beta_slope, capcal_slope,&
    7979                           albedo_slope,emiss_slope,qsurf_slope,&
    80                            iflat,ini_comslope_h
     80                           iflat,major_slope,ini_comslope_h
    8181      use time_phylmdz_mod, only: daysec,dtphys
    8282      USE comconst_mod, ONLY: rad,g,r,cpp,pi
    … …  
    272272     REAL,ALLOCATABLE  :: alph_locslope(:,:)
    273273     REAL,ALLOCATABLE  :: beta_locslope(:,:)   
    274 
     274     REAL,ALLOCATABLE  :: watersurf_density_timeseries(:,:,:,:)
     275     REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:,:)
     276     REAL,ALLOCATABLE  :: watersurf_density_phys_timeseries(:,:,:)
     277     REAL,ALLOCATABLE  :: watersurf_density_phys_ave(:,:)
     278     REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_timeseries(:,:,:,:)
     279     REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_ave(:,:,:)
    275280     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)
    276281     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)
    277282     REAL :: totmass_adsorbded
     283   real :: alpha_clap_h2o = -6143.7
     284   real :: beta_clap_h2o = 28.9074
    278285
    279286#ifdef CPP_STD
    … …  
    289296! Loop variable
    290297     LOGICAL :: bool_sublim
    291      INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop
     298     INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil
    292299     REAL :: beta = 3182.48
    293300     REAL :: alpha = 23.3494
    … …  
    406413              watercap,inertiesoil,nslope,tsurf_slope,           &
    407414              tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
    408               subslope_dist,albedo_slope,emiss_slope, TI_GCM_start,     &
     415              subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM_start,     &
    409416              qsurf_slope,watercap_slope)
    410417
    … …  
    551558     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    552559     allocate(delta_co2_adsorbded(ngrid))
    553 
     560     allocate(watersurf_density_timeseries(iim+1,jjm+1,nslope,timelen))
     561     allocate(watersoil_density_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
     562     allocate(watersurf_density_phys_timeseries(ngrid,nslope,timelen))
     563     allocate(watersurf_density_phys_ave(ngrid,nslope))
     564     allocate(watersoil_density_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
     565     allocate(watersoil_density_phys_PEM_ave(ngrid,nsoilmx_PEM,nslope))
    554566     print *, "Downloading data Y1..."
    555567
    556568     call read_data_GCM("data_GCM_Y1.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_yr1,timelen,min_co2_ice_slope_1,min_h2o_ice_slope_1,&   
    557                        nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope)
     569                       nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &     
     570                       watersurf_density_timeseries,watersoil_density_timeseries)
    558571
    559572! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value
    … …  
    569582
    570583     call read_data_GCM("data_GCM_Y2.nc", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM,timelen,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
    571                   nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope)
     584                  nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &     
     585                  watersurf_density_timeseries,watersoil_density_timeseries)
    572586
    573587     print *, "Downloading data Y2 done"
    … …  
    622636
    623637      allocate(ice_depth(ngrid,nslope))
     638      ice_depth(:,:) = 0.
    624639      allocate(TI_GCM_phys(ngrid,nsoilmx,nslope))
    625640
    … …  
    636651      deallocate(co2_ice_GCM_slope)
    637652      deallocate(TI_GCM)
    638 
    639   if(soil_pem) then
    640 
    641653      deallocate(tsurf_GCM_timeseries)
     654
     655
     656   if(soil_pem) then
    642657      call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
    643 
    644658      DO islope = 1,nslope
     659        DO t=1,timelen
     660          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersurf_density_timeseries(:,:,islope,t),watersurf_density_phys_timeseries(:,islope,t))
     661        ENDDO
    645662        DO l=1,nsoilmx
    646663           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
    … …  
    648665           DO t=1,timelen
    649666             CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
     667           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersoil_density_timeseries(:,:,l,islope,t),watersoil_density_phys_PEM_timeseries(:,l,islope,t))
    650668           ENDDO
     669
    651670        ENDDO
    652671        DO l=nsoilmx+1,nsoilmx_PEM
    653672           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
    654673           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
     674           DO t=1,timelen
     675           watersoil_density_phys_PEM_timeseries(:,l,islope,t) =  watersoil_density_phys_PEM_timeseries(:,nsoilmx,islope,t)
     676           ENDDO
    655677        ENDDO
    656678      ENDDO
    657 
     679       watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
     680      watersurf_density_phys_ave(:,:) = SUM(watersurf_density_phys_timeseries(:,:,:),3)/timelen
     681      deallocate(watersurf_density_timeseries)
     682      deallocate(watersurf_density_phys_timeseries)
     683      deallocate(watersoil_density_timeseries)
    658684      deallocate(tsoil_ave_yr1)
    659685      deallocate(tsoil_ave)
    … …  
    789815!---------------------------- Read the PEMstart ---------------------
    790816
    791       call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
     817      call pemetat0(.true.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
    792818      tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,&
    793       tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), co2_adsorbded_phys,delta_co2_adsorbded)
     819      tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM, co2_adsorbded_phys,delta_co2_adsorbded,&
     820      watersurf_density_phys_ave,watersoil_density_phys_PEM_ave)
    794821
    795822    if(soil_pem) then
    … …  
    848875     year_iter=0
    849876     print *, "Max_iter_pem", Max_iter_pem
     877     print *, 'year_iter_max',year_iter_max
    850878     do while (year_iter.LT.year_iter_max)
    851879
    … …  
    856884           global_ave_press_new=global_ave_press_new-g*cell_area(i)*tendencies_co2_ice_phys_slope(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
    857885       enddo
     886     enddo
     887       print *, 'Global average pressure old time step',global_ave_press_old
     888       print *, 'Global average pressure new time step',global_ave_press_new
     889
     890     do i=1,ngrid
    858891       if(soil_pem) then
    859892         global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
    860893       endif
    861894     enddo
     895       print *, 'Global average pressure old time step',global_ave_press_old
     896       print *, 'Global average pressure new time step',global_ave_press_new
    862897
    863898! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
    … …  
    9741009              flag_co2flow_mesh(ig) = 1.
    9751010            ENDIF ! co2ice > hmax
    976           ENDIF ! iflattsoil_lope
     1011          ENDIF ! iflat
    9771012        ENDDO !islope
    9781013       ENDDO !ig
    … …  
    10321067              emiss_slope(ig,islope) = emissiv         
    10331068            endif
    1034           elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2
     1069          elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice_phys_slope(ig,islope).gt.1e-5) )THEN !Put tsurf as tcond co2
    10351070            ave=0.
    10361071            do t=1,timelen
    1037               if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then
     1072              if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-3) then
    10381073                ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
    10391074              else
    … …  
    10561091      enddo
    10571092
     1093
    10581094    if(soil_pem) then
    10591095
    … …  
    10711107  do islope = 1,nslope
    10721108     TI_locslope(:,:) = TI_PEM(:,:,islope)
    1073      Tsoil_locslope(:,:) = tsoil_PEM(:,:,islope)
    1074      Tsurf_locslope(:) = tsurf_ave_phys(:,islope)
    10751109     alph_locslope(:,:) = alph_PEM(:,:,islope)
    10761110     beta_locslope(:,:) = beta_PEM(:,:,islope)
    1077      call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
    1078      call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
    10791111     do t = 1,timelen
    1080        tsoil_phys_PEM_timeseries(:,:,islope,t) =  tsoil_phys_PEM_timeseries(:,:,islope,t) + (Tsoil_locslope(:,:)- tsoil_PEM(:,:,islope))
    1081      enddo
    1082      TI_PEM(:,:,islope) = TI_locslope(:,:)
    1083      tsoil_PEM(:,:,islope) = Tsoil_locslope(:,:)
    1084      tsurf_ave_phys(:,islope)  =  Tsurf_locslope(:)
    1085      alph_PEM(:,:,islope) = alph_locslope(:,:)
    1086      beta_PEM(:,:,islope) = beta_locslope(:,:)
    1087   enddo
    1088 
    1089   print *, "Update of soil temperature done"
    1090 
    1091 ! Check Nan in the time series
    1092      do ig = 1,ngrid
    1093        do islope = 1,nslope
    1094          do iloop = 1,nsoilmx_PEM
    1095            if(isnan(tsoil_PEM(ig,iloop,islope))) then
    1096             write(*,*) "Problem : There is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope)
    1097            stop
    1098            endif
     1112       Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
     1113       Tsurf_locslope(:) =  tsurf_phys_GCM_timeseries(:,islope,t)
     1114
     1115       call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
     1116       call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
     1117
     1118
     1119       tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
     1120       do ig = 1,ngrid
     1121         do isoil = 1,nsoilmx_PEM
     1122          watersoil_density_phys_PEM_timeseries(ig,l,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
     1123          if(isnan(Tsoil_locslope(ig,isoil))) then
     1124             write(*,*)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil)
     1125          endif
    10991126         enddo
    11001127       enddo
    11011128     enddo
     1129     alph_PEM(:,:,islope) = alph_locslope(:,:)
     1130     beta_PEM(:,:,islope) = beta_locslope(:,:)
     1131  enddo
     1132
     1133 
     1134     tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
     1135     watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
     1136  print *, "Update of soil temperature done"
     1137
     1138  deallocate(TI_locslope)
     1139  deallocate(Tsoil_locslope)
     1140  deallocate(Tsurf_locslope)
     1141  deallocate(alph_locslope)
     1142  deallocate(beta_locslope)
    11021143
    11031144  write(*,*) "Compute ice table"
    11041145
    11051146! II_d.3 Update the ice table
    1106      call computeice_table(timelen,ngrid,nslope,nsoilmx,nsoilmx_PEM, tsoil_phys_PEM_timeseries,tsurf_phys_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
    1107                            ps_phys_timeseries, ice_depth)
     1147     call computeice_table(ngrid,nslope,nsoilmx_PEM,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
    11081148       
    11091149      print *, "Update soil propreties"
    … …  
    11531193      if (STOPPING_water) then
    11541194        print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
     1195      endif
     1196      if (year_iter.ge.year_iter_max) then
     1197        print *, "STOPPING because maximum number of iterations reached"
    11551198      endif
    11561199      if (STOPPING_1_water) then
    … …  
    12101253      ENDDO ! of DO ig=1,ngrid
    12111254
    1212 ! III_a.2 Tsurf and Tsoil update (for startfi)
     1255! III_a.2 Tsoil update (for startfi)
    12131256
    12141257   if(soil_pem) then
    … …  
    12191262   endif !soil_pem
    12201263
    1221       DO ig = 1,ngrid
    1222         tsurf(ig) = 0.
    1223         DO islope = 1,nslope
    1224           tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) &
    1225                      * subslope_dist(ig,islope)     
    1226         ENDDO
    1227       ENDDO ! of DO ig=1,ngrid
    12281264
    12291265#ifndef CPP_STD
    … …  
    12601296      ENDDO
    12611297    ENDDO
     1298
     1299
     1300      DO ig = 1,ngrid
     1301        tsurf(ig) = 0.
     1302        DO islope = 1,nslope
     1303          tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)*tsurf_slope(ig,islope))**4 &
     1304                     * subslope_dist(ig,islope)     
     1305        ENDDO
     1306        tsurf(ig) = tsurf(ig)**(0.25)/emis(ig)
     1307      ENDDO ! of DO ig=1,ngrid
     1308
    12621309#endif
    12631310

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

    r2835 r2849  
    11   SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
    2                       tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, m_co2_regolith_phys,deltam_co2_regolith_phys)
     2                      tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
     3                      global_ave_pressure, m_co2_regolith_phys,deltam_co2_regolith_phys, watersurf_ave,watersoil_ave)
    34   
    45   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
    56   use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
    67   use comsoil_h, only:  volcapa,inertiedat
    7    use ini_soil_mod, only:  ini_icetable
    88   use soil_evolution_mod, only: soil_pem,soil_pem_CN
    99   use adsorption_mod, only : regolith_co2adsorption
    … …  
    3232  real,intent(in) :: waterice(ngrid,nslope)
    3333  real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope)
    34 
     34  real, intent(in) :: watersurf_ave(ngrid,nslope)     ! surface water ice density, yearly averaged  (kg/m^3)
     35  real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)    ! surface water ice density, yearly averaged (kg/m^3)
     36  real, intent(in) :: global_ave_pressure                       ! mean average pressure on the planet
    3537! outputs
    3638
    … …  
    5961   real :: co2_ads_prev(ngrid)
    6062   real :: year_PEM_read
     63   real :: alpha_clap_h2o = -6143.7
     64   real :: beta_clap_h2o = 28.9074
    6165
    6266!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    221225        enddo
    222226     enddo
     227      do isoil = nsoil_GCM+1,nsoil_PEM
     228        do ig = 1,ngrid
     229        watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope)
     230        enddo
     231      enddo
     232
    223233   
    224234ENDDO
    … …  
    230240!3. Ice Table
    231241
    232    call get_field("ice_table",ice_table,found)
    233    if(.not.found) then
    234       write(*,*)'PEM settings: failed loading <Ice Table>'
    235       write(*,*)'will reconstruct the values of ice table'
    236 
    237       call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:),  tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
    238 
    239    else
    240    ! update ice table
    241      call computeice_table(timelen,ngrid,nslope,nsoil_PEM,tsoil_inst,tsurf_inst,q_co2,q_h2o, ps_inst, ice_table)
    242 
    243    endif
     242
     243
     244     call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
     245     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)
     246
    244247
    245248print *,'PEMETAT0: ICE TABLE  DONE'
    … …  
    259262   deltam_co2_regolith_phys(:) = 0.
    260263   exit
     264   else
     265   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
     266
    261267   endif
    262268ENDDO
    263269
    264 if (found) then
    265    DO iloop = 1,nsoil_GCM
    266       tsoil_tmp_yr1(:,iloop,:) = tsoil_PEM_yr1(:,iloop,:)
    267 
    268    ENDDO 
    269 
    270    call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
    271 
    272 
    273 endif
     270
    274271print *,'PEMETAT0: CO2 adsorption done  '
    275272
    … …  
    364361!b) Soil temperature
    365362    do islope = 1,nslope
    366 
    367      write(*,*) "islope=",islope
    368363     do ig = 1,ngrid
    369364        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
    … …  
    381376            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
    382377      enddo
    383 !     write(*,*) "ig, islope, T=", ig,islope,tsoil_PEM(ig,:,islope)
    384378     enddo
    385379   
    … …  
    389383        enddo
    390384     enddo
     385      do isoil = nsoil_GCM+1,nsoil_PEM
     386        do ig = 1,ngrid
     387        watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope)
     388        enddo
     389      enddo
    391390enddo
    392391print *,'PEMETAT0: TSOIL DONE  '
    393392!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    394393!c) Ice table   
    395       do islope = 1,nslope
    396       call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:), & 
    397        tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
    398       enddo
     394     call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
     395     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)
     396
     397
    399398!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    400399

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

    r2842 r2849  
    33!
    44SUBROUTINE read_data_GCM(fichnom,min_h2o_ice_s,min_co2_ice_s,iim_input,jjm_input,nlayer,vmr_co2_gcm,ps_GCM,timelen, &
    5              min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope)
     5             min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope, &
     6             watersurf_density,watersoil_density)
    67
    78      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
    … …  
    3940  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,timelen)
    4041  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,timelen)
    41   REAL, ALLOCATABLE ::  q1_co2_GCM(:,:,:)
    4242  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,timelen)
    4343
    … …  
    4848  REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file
    4949  REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
     50  REAL , INTENT(OUT) ::  watersurf_density(iim_input+1,jjm_input+1,nslope,timelen) ! Soil temperature of the concatenated file
     51  REAL , INTENT(OUT) ::  watersoil_density(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
    5052
    5153  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Thermal Inertia  of the concatenated file
    … …  
    7981
    8082      allocate(co2_ice_s(iim+1,jjm+1,timelen))
    81       allocate(q1_co2_GCM(iim+1,jjm+1,timelen))
    8283      allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
    8384      allocate(h2o_ice_s(iim+1,jjm+1,timelen))
    … …  
    161162
    162163     print *, "Downloading data for inertiesoil_slope done"
     164
     165     print *, "Downloading data for watersoil_density ..."
     166
     167DO islope=1,nslope
     168  write(num,fmt='(i2.2)') islope
     169  call get_var4("Waterdensity_soil_slope"//num,watersoil_density(:,:,:,islope,:))
     170ENDDO
     171
     172     print *, "Downloading data for  watersoil_density  done"
     173
     174     print *, "Downloading data for  watersurf_density  ..."
     175
     176DO islope=1,nslope
     177  write(num,fmt='(i2.2)') islope
     178  call get_var3("Waterdensity_surface"//num,watersurf_density(:,:,islope,:))
     179ENDDO
     180
     181     print *, "Downloading data for  watersurf_density  done"
    163182
    164183  endif !soil_pem
    … …  
    234253  ENDDO
    235254
     255
     256      deallocate(co2_ice_s)
     257      deallocate(h2o_ice_s_slope)
     258      deallocate(h2o_ice_s)
     259
     260
    236261  CONTAINS
    237262

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

    r2842 r2849  
    88
    99      USE temps_mod_evol, ONLY: year_bp_ini, year_PEM
    10 #ifndef CPP_STD
     10#ifndef CPP_STD     
     11      USE comconst_mod, ONLY: pi
    1112      USE planete_h, ONLY: e_elips, obliquit, timeperi
    1213#else
    1314      use planete_mod, only: e_elips, obliquit, timeperi
     15      USE comcstfi_mod, only: pi
     16
    1417#endif
    15       USE comcstfi_mod, only: pi
     18
    1619      USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
    1720                                end_lask_param_mod,last_ilask

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

    r2842 r2849  
    2626 REAL ::  L =  51058.
    2727 REAL ::  inertie_thresold = 800. ! look for ice
    28  REAL ::  inertie_co2glaciers = 2120 ! Mellon et al. 2000
    2928 REAL ::  inertie_averaged = 250 ! Mellon et al. 2000
    30  REAL ::  ice_inertia = 2120  ! Inertia of ice
    31  REAL ::  surfaceice_inertia = 800  ! Inertia of ice
     29 REAL ::  ice_inertia = 1200  ! Inertia of ice
    3230 REAL ::  P610 = 610.0
    3331 REAL ::  m_h2o = 18.01528E-3
    … …  
    4745! 0. Initialisation
    4846
    49   do ig = 1,ngrid
    50     do islope = 1,nslope
    51      if((abs(tend_h2oglaciers(ig,islope)).lt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then
    52         ispermanent_h2oglaciers(ig,islope) = 1
    53      else
    54         ispermanent_h2oglaciers(ig,islope) = 0
    55      endif
     47!  do ig = 1,ngrid
     48!    do islope = 1,nslope
     49!     if((abs(tend_h2oglaciers(ig,islope)).lt.1e-5).and.(abs(waterice(ig,islope)).gt.(1.e-4))) then
     50!        ispermanent_h2oglaciers(ig,islope) = 1
     51!     else
     52!        ispermanent_h2oglaciers(ig,islope) = 0
     53!     endif
     54!
     55!     if((abs(tend_co2glaciers(ig,islope)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.(1.e-4))) then
     56!        ispermanent_co2glaciers(ig,islope) = 1
     57!     else
     58!        ispermanent_co2glaciers(ig,islope) = 0
     59!     endif
     60!    enddo
     61!  enddo
    5662
    57      if((abs(tend_co2glaciers(ig,islope)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
    58         ispermanent_co2glaciers(ig,islope) = 1
    59      else
    60         ispermanent_co2glaciers(ig,islope) = 0
    61      endif
    62     enddo
    6363
    64   enddo
    65 
    66  ispermanent_co2glaciers(:,:) = 0
    67  ispermanent_h2oglaciers(:,:) = 0
    6864
    6965!  1.Ice TI feedback
    7066
    71     do islope = 1,nslope
    72      call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope),   TI_PEM(:,:,islope))
    73     enddo
     67!    do islope = 1,nslope
     68!     call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope),   TI_PEM(:,:,islope))
     69!    enddo
    7470
    7571! 2. Modification of the regolith thermal inertia.
Note: See TracChangeset for help on using the changeset viewer.