Ignore:
Timestamp:
Oct 6, 2023, 5:32:11 PM (14 months ago)
Author:
jbclement
Message:

PEM:
Big cleaning/improvements of the PEM:

  • Conversion of "abort_pem.F" and "soil_settings_PEM.F" into Fortran 90;
  • Transformation of every PEM subroutines into module;
  • Rewriting of many subroutines with modern Fortran syntax;
  • Correction of a bug in "pem.F90" when calling 'recomp_tend_co2_slope'. The arguments were given in disorder and emissivity was missing;
  • Update of "launch_pem.sh" in deftank.

JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3039 r3076  
    1 SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
    2                       tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
    3                       global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
    4                       m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
    5    
    6   use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
    7   use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
    8   use comsoil_h,                  only: volcapa,inertiedat
    9   use adsorption_mod,             only: regolith_adsorption, adsorption_pem
    10   use ice_table_mod,              only: computeice_table_equilibrium
    11   use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
    12   use soil_thermalproperties_mod, only: update_soil_thermalproperties
    13   use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
     1MODULE pemetat0_mod
     2
     3implicit none
     4
     5!=======================================================================
     6contains
     7!=======================================================================
     8
     9SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness,   &
     10                    tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
     11                    global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,                &
     12                    m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
     13
     14use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
     15use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
     16use comsoil_h,                  only: volcapa,inertiedat
     17use adsorption_mod,             only: regolith_adsorption, adsorption_pem
     18use ice_table_mod,              only: computeice_table_equilibrium
     19use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
     20use soil_thermalproperties_mod, only: update_soil_thermalproperties
     21use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
     22use abort_pem_mod,              only: abort_pem
     23use soil_pem_compute_mod,       only: soil_pem_compute
     24use soil_pem_ini_mod,           only: soil_pem_ini
    1425
    1526#ifndef CPP_STD
    … …  
    2031#endif
    2132
    22   implicit none
    23 
    24   include "callkeys.h"
    25 
    26   character(len=*), intent(in) :: filename                    ! name of the startfi_PEM.nc
    27   integer,intent(in) :: ngrid                                 ! # of physical grid points
    28   integer,intent(in) :: nsoil_GCM                             ! # of vertical grid points in the GCM
    29   integer,intent(in) :: nsoil_PEM                             ! # of vertical grid points in the PEM
    30   integer,intent(in) :: nslope                                ! # of sub-grid slopes
    31   integer,intent(in) :: timelen                               ! # time samples
    32   real, intent(in) :: tsurf_ave_yr1(ngrid,nslope)             ! surface temperature at the first year of GCM call [K]
    33   real,intent(in) :: tsurf_ave_yr2(ngrid,nslope)              ! surface temperature at the second  year of GCM call [K]
    34   real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
    35   real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
    36   real,intent(in) :: ps_inst(ngrid,timelen)                   ! surface pressure [Pa]
    37   real,intent(in) :: timestep                                 ! time step [s]
    38   real,intent(in) :: tend_h2oglaciers(ngrid,nslope)           ! tendencies on h2o glaciers
    39   real,intent(in) :: tend_co2glaciers(ngrid,nslope)           ! tendencies on co2 glaciers
    40   real,intent(in) :: co2ice(ngrid,nslope)                     ! co2 ice amount [kg/m^2]
    41   real,intent(in) :: waterice(ngrid,nslope)                   ! water ice amount [kg/m^2]
    42   real, intent(in) :: watersurf_ave(ngrid,nslope)             ! surface water ice density, yearly averaged  (kg/m^3)
    43   real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)! surface water ice density, yearly averaged (kg/m^3)
    44   real, intent(in) :: global_ave_pressure                     ! mean average pressure on the planet [Pa]
     33implicit none
     34
     35include "callkeys.h"
     36
     37character(len=*),               intent(in) :: filename            ! name of the startfi_PEM.nc
     38integer,                        intent(in) :: ngrid               ! # of physical grid points
     39integer,                        intent(in) :: nsoil_GCM           ! # of vertical grid points in the GCM
     40integer,                        intent(in) :: nsoil_PEM           ! # of vertical grid points in the PEM
     41integer,                        intent(in) :: nslope              ! # of sub-grid slopes
     42integer,                        intent(in) :: timelen             ! # time samples
     43real,                           intent(in) :: timestep            ! time step [s]
     44real,                           intent(in) :: global_ave_pressure ! mean average pressure on the planet [Pa]
     45real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr1       ! surface temperature at the first year of GCM call [K]
     46real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr2       ! surface temperature at the second  year of GCM call [K]
     47real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
     48real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
     49real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
     50real, dimension(ngrid,nslope),  intent(in) :: tend_h2oglaciers    ! tendencies on h2o glaciers
     51real, dimension(ngrid,nslope),  intent(in) :: tend_co2glaciers    ! tendencies on co2 glaciers
     52real, dimension(ngrid,nslope),  intent(in) :: co2ice              ! co2 ice amount [kg/m^2]
     53real, dimension(ngrid,nslope),  intent(in) :: waterice            ! water ice amount [kg/m^2]
     54real, dimension(ngrid,nslope),  intent(in) :: watersurf_ave       ! surface water ice density, yearly averaged  (kg/m^3)
    4555! outputs
    46 
    47   real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope)                ! soil (mid-layer) thermal inertia in the PEM grid [SI]
    48   real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)             ! soil (mid-layer) temperature [K]
    49   real,intent(inout) :: ice_table(ngrid,nslope)                       ! Ice table depth [m]
    50   real,intent(inout) :: ice_table_thickness(ngrid,nslope)             ! Ice table thickness [m]
    51   real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen)    ! instantaneous soil (mid-layer) temperature [K]
    52   real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed [kg/m^2]
    53   real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
    54   real,intent(inout) :: m_h2o_regolith_phys(ngrid,nsoil_PEM,nslope)  ! mass of h2o adsorbed [kg/m^2]
    55   real,intent(out) :: deltam_h2o_regolith_phys(ngrid)                 ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
    56   real,intent(inout) :: water_reservoir(ngrid)                        ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
     56real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
     57real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
     58real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
     59real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
     60real, dimension(ngrid,nslope),                   intent(inout) :: ice_table           ! Ice table depth [m]
     61real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
     62real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
     63real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
     64real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
     65real, dimension(ngrid),                          intent(inout) :: water_reservoir     ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
     66real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_ave       ! surface water ice density, yearly averaged (kg/m^3)
    5767! local
    58    real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)                     ! soil temperature saved in the start [K]
    59    real :: TI_startPEM(ngrid,nsoil_PEM,nslope)                        ! soil thermal inertia saved in the start [SI]
    60    LOGICAL :: found                                                   ! check if variables are found in the start
    61    LOGICAL :: found2                                                  ! check if variables are found in the start
    62    integer :: iloop,ig,islope,it,isoil                                ! index for loops
    63    real :: kcond                                                      ! Thermal conductivity, intermediate variable [SI]
    64    real :: delta                                                      ! Depth of the interface regolith-breccia, breccia -bedrock [m]
    65    CHARACTER*2 :: num                                                 ! intermediate string to read PEM start sloped variables
    66    real :: tsoil_saved(ngrid,nsoil_PEM)                               ! saved soil temperature [K]
    67    real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope)                      ! intermediate soil temperature during yr1[K]
    68    real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)                      ! intermediate soil temperature during yr 2 [K]
    69    real :: alph_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computation []
    70    real :: beta_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computatio []                               
    71    LOGICAL :: startpem_file                                           ! boolean to check if we read the startfile or not
    72 #ifdef CPP_STD   
    73    logical :: watercaptag(ngrid)
    74 
    75    watercaptag(:)=.false.
     68real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM               ! soil temperature saved in the start [K]
     69real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                  ! soil thermal inertia saved in the start [SI]
     70logical                                 :: found                        ! check if variables are found in the start
     71logical                                 :: found2                       ! check if variables are found in the start
     72integer                                 :: iloop, ig, islope, it, isoil ! index for loops
     73real                                    :: kcond                        ! Thermal conductivity, intermediate variable [SI]
     74real                                    :: delta                        ! Depth of the interface regolith-breccia, breccia -bedrock [m]
     75character(2)                            :: num                          ! intermediate string to read PEM start sloped variables
     76real, dimension(ngrid,nsoil_PEM)        :: tsoil_saved                  ! saved soil temperature [K]
     77real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr1[K]
     78real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                ! intermediate soil temperature during yr 2 [K]
     79real, dimension(ngrid,nsoil_PEM - 1)    :: alph_tmp                     ! Intermediate for tsoil computation []
     80real, dimension(ngrid,nsoil_PEM - 1)    :: beta_tmp                     ! Intermediate for tsoil computatio []
     81logical                                 :: startpem_file                ! boolean to check if we read the startfile or not
     82
     83#ifdef CPP_STD
     84   logical, dimension(ngrid) :: watercaptag
     85   watercaptag(:) = .false.
    7686#endif
    7787
    … …  
    8191!!!
    8292!!! Order: 0. Previous year of the PEM run
    83 !!!        1. Thermal Inertia 
     93!!!        1. Thermal Inertia
    8494!!!        2. Soil Temperature
    8595!!!        3. Ice table
    … …  
    99109
    100110!1. Run
    101 
    102111if (startpem_file) then
    103    ! open pem initial state file:
    104    call open_startphy(filename)
     112    ! open pem initial state file:
     113    call open_startphy(filename)
    105114
    106115    if (soil_pem) then
    107  
     116
    108117!1. Thermal Inertia
    109118! a. General case
    110 
    111 DO islope=1,nslope
    112   write(num,fmt='(i2.2)') islope
    113   call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
    114    if(.not.found) then
    115       write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
    116       write(*,*)'will reconstruct the values of TI_PEM'
    117 
    118       do ig = 1,ngrid
    119           if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then
    120 !!! transition
    121              delta = depth_breccia
    122              TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
    123             (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
    124                       ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
    125             do iloop=index_breccia+2,index_bedrock
    126               TI_PEM(ig,iloop,islope) = TI_breccia
    127             enddo     
    128           else ! we keep the high ti values
    129             do iloop=index_breccia+1,index_bedrock
    130               TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
    131             enddo
    132           endif ! TI PEM and breccia comparison
    133 !! transition
    134           delta = depth_bedrock
    135           TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
    136                (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
    137                ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
    138           do iloop=index_bedrock+2,nsoil_PEM
    139             TI_PEM(ig,iloop,islope) = TI_bedrock
    140           enddo   
    141       enddo
    142    else ! found
    143      do iloop = nsoil_GCM+1,nsoil_PEM
    144        TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
    145      enddo
    146    endif ! not found
    147 ENDDO ! islope
    148 
    149 write(*,*) 'PEMETAT0: THERMAL INERTIA DONE'
     119    do islope = 1,nslope
     120        write(num,'(i2.2)') islope
     121        call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
     122        if (.not. found) then
     123            write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
     124            write(*,*)'will reconstruct the values of TI_PEM'
     125
     126            do ig = 1,ngrid
     127                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
     128                    !!! transition
     129                    delta = depth_breccia
     130                    TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
     131                                                        (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
     132                                                        ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
     133                    do iloop=index_breccia+2,index_bedrock
     134                        TI_PEM(ig,iloop,islope) = TI_breccia
     135                    enddo
     136                else ! we keep the high ti values
     137                    do iloop = index_breccia + 1,index_bedrock
     138                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
     139                    enddo
     140                endif ! TI PEM and breccia comparison
     141                !! transition
     142                delta = depth_bedrock
     143                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
     144                                                        (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
     145                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
     146                do iloop = index_bedrock + 2,nsoil_PEM
     147                    TI_PEM(ig,iloop,islope) = TI_bedrock
     148                enddo
     149            enddo
     150        else ! found
     151            do iloop = nsoil_GCM + 1,nsoil_PEM
     152                TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
     153            enddo
     154        endif ! not found
     155    enddo ! islope
     156
     157    write(*,*) 'PEMETAT0: THERMAL INERTIA doNE'
    150158
    151159! b. Special case for inertiedat, inertiedat_PEM
    … …  
    154162 do iloop = 1,nsoil_GCM
    155163   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
    156  enddo 
     164 enddo
    157165!!! zone de transition
    158166delta = depth_breccia
    159167do ig = 1,ngrid
    160 if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then 
     168if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
    161169inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
    162170            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
    163171                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
    164172
    165  do iloop = index_breccia+2,index_bedrock 
     173 do iloop = index_breccia+2,index_bedrock
    166174       inertiedat_PEM(ig,iloop) = TI_breccia
    167175  enddo
    … …  
    170178   do iloop=index_breccia+1,index_bedrock
    171179      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
    172    enddo 
     180   enddo
    173181endif ! comparison ti breccia
    174182enddo!ig
    … …  
    189197endif ! not found
    190198
    191 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     199!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    192200!2. Soil Temperature
    193 DO islope=1,nslope
     201do islope=1,nslope
    194202  write(num,fmt='(i2.2)') islope
    195203   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
    … …  
    199207!      do ig = 1,ngrid
    200208!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
    201 !        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))   
     209!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
    202210!       do iloop=index_breccia+2,index_bedrock
    203211!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
    … …  
    205213!      enddo
    206214!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
    207 !        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1)) 
     215!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
    208216!
    209217!      do iloop=index_bedrock+2,nsoil_PEM
    … …  
    215223
    216224     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    217      call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    218 
    219    
     225     call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     226
     227
    220228   else
    221229! predictor corrector: restart from year 1 of the GCM and build the evolution of
    … …  
    223231
    224232    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
    225     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
    226     call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
    227     tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
    228     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
     233    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
     234    call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
     235    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
     236    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
    229237
    230238
    … …  
    244252        enddo
    245253      enddo
    246 ENDDO ! islope
    247 
    248 write(*,*) 'PEMETAT0: SOIL TEMP  DONE'
    249 
    250 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     254enddo ! islope
     255
     256write(*,*) 'PEMETAT0: SOIL TEMP  doNE'
     257
     258!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    251259!3. Ice Table
    252260  call get_field("ice_table",ice_table,found)
    … …  
    261269   endif
    262270
    263 write(*,*) 'PEMETAT0: ICE TABLE  DONE'
    264 
    265 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     271write(*,*) 'PEMETAT0: ICE TABLE  doNE'
     272
     273!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    266274!4. CO2 & H2O Adsorption
    267275 if(adsorption_pem) then
    268   DO islope=1,nslope
     276  do islope=1,nslope
    269277   write(num,fmt='(i2.2)') islope
    270278   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
    … …  
    273281       exit
    274282    endif
    275    
    276   ENDDO
    277 
    278   DO islope=1,nslope
     283
     284  enddo
     285
     286  do islope=1,nslope
    279287   write(num,fmt='(i2.2)') islope
    280288   call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
    … …  
    283291      exit
    284292    endif
    285    
    286   ENDDO
     293
     294  enddo
    287295
    288296    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
    … …  
    293301    endif
    294302    if((.not.found2)) then
    295        deltam_h2o_regolith_phys(:) = 0. 
     303       deltam_h2o_regolith_phys(:) = 0.
    296304    endif
    297 write(*,*) 'PEMETAT0: CO2 & H2O adsorption DONE'
     305write(*,*) 'PEMETAT0: CO2 & H2O adsorption doNE'
    298306    endif
    299307endif ! soil_pem
    300308
    301 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     309!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    302310!. 5 water reservoir
    303 #ifndef CPP_STD   
     311#ifndef CPP_STD
    304312   call get_field("water_reservoir",water_reservoir,found)
    305313   if(.not.found) then
    … …  
    326334      do ig = 1,ngrid
    327335
    328           if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then 
     336          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then
    329337!!! transition
    330338             delta = depth_breccia
    … …  
    335343          do iloop=index_breccia+2,index_bedrock
    336344            TI_PEM(ig,iloop,islope) = TI_breccia
    337          enddo     
     345         enddo
    338346         else ! we keep the high ti values
    339347           do iloop=index_breccia+1,index_bedrock
    340348                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
    341            enddo 
     349           enddo
    342350         endif
    343351
    … …  
    349357          do iloop=index_bedrock+2,nsoil_PEM
    350358            TI_PEM(ig,iloop,islope) = TI_bedrock
    351          enddo   
     359         enddo
    352360      enddo
    353361enddo
    … …  
    357365 enddo
    358366!!! zone de transition
    359 delta = depth_breccia 
     367delta = depth_breccia
    360368do ig = 1,ngrid
    361       if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then 
     369      if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
    362370inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
    363371            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
    … …  
    365373
    366374
    367  do iloop = index_breccia+2,index_bedrock 
     375 do iloop = index_breccia+2,index_bedrock
    368376
    369377       inertiedat_PEM(ig,iloop) = TI_breccia
    370    
     378
    371379 enddo
    372 else 
    373    do iloop = index_breccia+1,index_bedrock 
     380else
     381   do iloop = index_breccia+1,index_bedrock
    374382       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
    375383    enddo
    … …  
    394402 enddo
    395403
    396 write(*,*) 'PEMETAT0: TI DONE'
    397 
    398 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    399 !b) Soil temperature 
     404write(*,*) 'PEMETAT0: TI doNE'
     405
     406!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     407!b) Soil temperature
    400408    do islope = 1,nslope
    401409!     do ig = 1,ngrid
    … …  
    414422!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
    415423!      enddo
    416      
     424
    417425!       enddo
    418426
    419427      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    420       call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    421 
    422    
     428      call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     429
     430
    423431    do it = 1,timelen
    424432        do isoil = nsoil_GCM+1,nsoil_PEM
    … …  
    433441        enddo
    434442enddo !islope
    435 write(*,*) 'PEMETAT0: TSOIL DONE'
    436 
    437 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     443write(*,*) 'PEMETAT0: TSOIL doNE'
     444
     445!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    438446!c) Ice table
    439447       call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
    … …  
    442450         call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    443451       enddo
    444        write(*,*) 'PEMETAT0: Ice table DONE'
    445 
    446 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     452       write(*,*) 'PEMETAT0: Ice table doNE'
     453
     454!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    447455!d) Regolith adsorbed
    448  if(adsorption_pem) then
     456if (adsorption_pem) then
    449457    m_co2_regolith_phys(:,:,:) = 0.
    450458    m_h2o_regolith_phys(:,:,:) = 0.
    451459
    452460    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
    453                                 m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
    454    
     461                             m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
     462
    455463    deltam_co2_regolith_phys(:) = 0.
    456     deltam_h2o_regolith_phys(:) = 0. 
    457  endif
    458 
    459 write(*,*) 'PEMETAT0: CO2 adsorption DONE'
     464    deltam_h2o_regolith_phys(:) = 0.
     465endif
     466
     467write(*,*) 'PEMETAT0: CO2 adsorption doNE'
    460468endif !soil_pem
    461469
    462 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     470!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    463471!. e) water reservoir
    464 #ifndef CPP_STD   
    465       do ig=1,ngrid
    466         if(watercaptag(ig)) then
    467            water_reservoir(ig)=water_reservoir_nom
     472#ifndef CPP_STD
     473    do ig = 1,ngrid
     474        if (watercaptag(ig)) then
     475            water_reservoir(ig) = water_reservoir_nom
    468476        else
    469            water_reservoir(ig)=0.
     477            water_reservoir(ig) = 0.
    470478        endif
    471       enddo
     479    enddo
    472480#endif
    473481
    474482endif ! of if (startphy_file)
    475483
    476 if(soil_pem) then
    477 !! Sanity check
    478     DO ig = 1,ngrid
    479       DO islope = 1,nslope
    480         DO iloop = 1,nsoil_PEM
    481           if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
    482         ENDDO
    483       ENDDO
    484      ENDDO
    485 endif!soil_pem
    486 
    487 END SUBROUTINE
     484if (soil_pem) then ! Sanity check
     485    do ig = 1,ngrid
     486        do islope = 1,nslope
     487            do iloop = 1,nsoil_PEM
     488                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
     489            enddo
     490        enddo
     491    enddo
     492endif !soil_pem
     493
     494END SUBROUTINE pemetat0
     495
     496END MODULE pemetat0_mod
     497
Note: See TracChangeset for help on using the changeset viewer.