Changeset 2835


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

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
3 added
4 deleted
23 edited

Legend:

Unmodified
Added
Removed

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

    r2794 r2835  
    1616 REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)
    1717 REAL,INTENT(INOUT) :: delta_mreg(ngrid)
    18 
    19 
    2018 
    2119! Constants:
    … …  
    7169enddo
    7270
    73 
    7471! 2. Check the exchange between the atmosphere and the regolith
    7572

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

    r2794 r2835  
    33  contains
    44
    5 
    6 
    75  subroutine regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, theta_h2o_adsorbded,m_h2o_adsorbed)
    86
    … …  
    119
    1210 implicit none
    13 
    1411
    1512 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! size dimension
    … …  
    2017 REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)   ! Soil temperature (K)
    2118
    22 
    2319! output
    2420 REAL,INTENT(OUT) ::  m_h2o_adsorbed(ngrid,nsoil_PEM,nslope)     ! Density of h2o adsorbed (kg/m^3)
    2521 REAL,INTENT(OUT) :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules
    26 
    2722
    2823! constant
    … …  
    4035 real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
    4136
    42 
    43 
    44 
    45 
    4637! local variable
    4738 real :: A,B                                                   ! Used to compute the mean mass above the surface
    … …  
    5445 real, allocatable :: pvapor_avg(:)                            ! yearly average vapor pressure above the surface
    5546
    56 
    57 
    5847! 0. Some initializations
    5948
    60 
    61    
    6249     allocate(mass_mean(ngrid,timelen))
    6350     allocate(zplev_mean(ngrid,timelen))
    … …  
    7764     enddo
    7865
    79 
    8066! b. pressure level
    8167     do it = 1,timelen
    … …  
    9278     deallocate(zplev_mean)
    9379     deallocate(mass_mean)
    94 
    95 
    96 
    97 
    9880
    9981! 2. we compute the mass of co2 adsorded in each layer of the meshes 
    … …  
    11092        theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1+K*p_sat))**mu
    11193        m_h2o_adsorbed(ig,iloop,islope) =as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith
    112       endif
    113  
     94     endif
    11495   enddo
    11596  enddo
    … …  
    11899 RETURN
    119100  end subroutine
    120 
    121101
    122102  SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_completesoil,delta_mreg)
    … …  
    140120 REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)   ! Density of co2 adsorbed (kg/m^3)
    141121 REAL,INTENT(INOUT) :: delta_mreg(ngrid)                            ! Difference density of co2 adsorbed (kg/m^3)
    142 
    143 
    144122 
    145123! Constants:
    … …  
    170148 real, allocatable :: pco2_avg(:)                              ! yearly averaged
    171149
    172 
    173 
    174150! 0. Some initializations
    175 
    176 
    177151
    178152     allocate(mass_mean(ngrid,timelen))
    … …  
    188162     delta_mreg(:) = 0.
    189163
    190 
    191 
    192164!0.1 Look at perenial ice
    193165  do ig = 1,ngrid
    … …  
    207179   enddo
    208180
    209 
    210181!   0.2  Compute the partial pressure of CO2
    211182!a. the molecular mass into the column
    … …  
    213184       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
    214185     enddo
    215 
    216186
    217187! b. pressure level
    … …  
    226196     pco2_avg(:) = sum(pco2(:,:),2)/timelen
    227197
    228  
    229198     deallocate(zplev_mean)
    230199     deallocate(mass_mean)
    … …  
    234203! 1. Compute the fraction of the pores occupied by H2O
    235204 call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,theta_h2o_adsorbed, m_h2o_adsorbed)
    236 
    237 
    238205
    239206! 2.  we compute the mass of co2 adsorded in each layer of the meshes 
    … …  
    256223 enddo
    257224enddo
    258 
    259225
    260226! 2. Check the exchange between the atmosphere and the regolith
    … …  
    278244   m_co2_completesoil(:,:,:) = dm_co2_regolith_slope(:,:,:)
    279245
    280 
    281246!=======================================================================
    282247      RETURN

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

    r2779 r2835  
    4747  ENDDO
    4848
    49 
    5049!  We reorganise the difference on the physical grid
    5150  tendencies_h2o_ice_phys(1)=tendencies_h2o_ice(1,1)
    … …  
    6160  tendencies_h2o_ice_phys(ig0) = tendencies_h2o_ice(1,jjm_input+1)
    6261
    63 
    6462END SUBROUTINE compute_tendencies
    6563

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

    r2794 r2835  
    66
    77      IMPLICIT NONE
    8 
    98
    109!=======================================================================
    … …  
    3433!=======================================================================
    3534
    36 
    3735!  We compute the difference
    38 !  tendencies_h2o_ice(:,:,:)=min_h2o_ice_Y2(:,:,:)-min_h2o_ice_Y1(:,:,:)
    3936
    4037  DO j=1,jjm_input+1
    … …  
    4643  ENDDO
    4744
    48      print *, "jjm_input+1", jjm_input+1
    49      print *, "iim_input+1", iim_input+1
    50      print *, "nslope+1", nslope+1
    51 
    5245!  If the difference is too small; there is no evolution
    5346  DO j=1,jjm_input+1
    5447    DO i = 1, iim_input
    5548       DO islope = 1, nslope
    56 !         print *, "tendencies_h2o_ice(i,j,islope)LAAA", tendencies_h2o_ice(i,j,islope)
    5749         if(abs(tendencies_h2o_ice(i,j,islope)).LT.1.0E-10) then
    5850            tendencies_h2o_ice(i,j,islope)=0.
    5951         endif
    60 !         print *, "tendencies_h2o_ice(i,j,islope)HERE", tendencies_h2o_ice(i,j,islope)
    6152       enddo
    6253    ENDDO
    6354  ENDDO
    6455
    65 
    66 !  We reorganise the difference on the physical grid
    67   tendencies_h2o_ice_phys(1,:)=tendencies_h2o_ice(1,1,:)
    68 
    69   ig0 = 2
    70   DO j=2,jjm_input
    71     DO i = 1, iim_input
    72        tendencies_h2o_ice_phys(ig0,:)  =tendencies_h2o_ice(i,j,:)
    73        ig0= ig0 + 1
    74     ENDDO
     56  DO islope = 1,nslope
     57    CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tendencies_h2o_ice(:,:,islope),tendencies_h2o_ice_phys(:,islope))
    7558  ENDDO
    76 
    77   tendencies_h2o_ice_phys(ig0,:) = tendencies_h2o_ice(1,jjm_input+1,:)
    78 
    79 !  print *, "tendencies_h2o_ice_physze", tendencies_h2o_ice_phys(:,:)
    80 
    8159
    8260END SUBROUTINE compute_tendencies_slope
    8361
    84 
    85 
    86 
    87 

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

    r2794 r2835  
    44    USE comsoil_h_PEM, only: layer_PEM
    55    implicit none
    6 
    76
    87    integer,intent(in) :: timelen,ngrid,nslope,nsoil_PEM,nsoil_GCM
    … …  
    1413    real,intent(out) :: ice_table(ngrid,nslope)                 ! ice table [m]
    1514
    16 
    1715    real :: m_h2o = 18.01528E-3
    1816    real :: m_co2 = 44.01E-3 
    … …  
    2119    real :: alpha = -6143.7
    2220    real :: beta = 29.9074
    23 
    2421
    2522    integer ig, islope,isoil,it
    … …  
    3734    real,allocatable :: diff_rho(:,:,:)                    ! difference of vapor content
    3835
    39 
    40 
    4136     allocate(rhovapor(ngrid,nslope,timelen))
    4237     allocate(rhovapor_avg(ngrid,nslope))
    … …  
    4540     allocate(mass_mean(ngrid,timelen))
    4641     allocate(zplev_mean(ngrid,timelen))
    47 
    48 
    49 
    5042
    5143! 0. Some initializations
    … …  
    6052       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
    6153     enddo
    62 
    6354
    6455! b. pressure level
    … …  
    9182     deallocate(pvapor)
    9283
    93 
    94 
    95 
    9684!    1.3 Density at the surface yearly averaged
    9785     rhovapor_avg(:,:) = SUM(rhovapor(:,:,:),3)/timelen
    9886
    99 
    10087     deallocate(rhovapor)
    101 
    102 
    103 
    10488
    10589! 2. Compute rho soil vapor
    … …  
    118102     enddo
    119103
    120 
    121 
    122104    rho_soil_avg(:,:,:) = SUM( rho_soil(:,:,:,:),4)/timelen
    123105    deallocate(rho_soil)
    124 
    125 
    126106
    127107! 3. Computing ice table
    … …  
    132112
    133113         do isoil = 1,nsoil_PEM
    134 
    135 
    136              diff_rho(:,:,isoil) = rhovapor_avg(:,:) - rho_soil_avg(:,:,isoil)
    137 !             write(*,*) 'diff =',ig,islope,isoil,diff_rho(ig,islope,isoil),rhovapor_avg(ig,islope) ,rho_soil_avg(ig,islope,isoil)
     114           diff_rho(:,:,isoil) = rhovapor_avg(:,:) - rho_soil_avg(:,:,isoil)
     115!          write(*,*) 'diff =',ig,islope,isoil,diff_rho(ig,islope,isoil),rhovapor_avg(ig,islope) ,rho_soil_avg(ig,islope,isoil)
    138116         enddo
    139117   
    140 
    141118  deallocate(rhovapor_avg)
    142       deallocate(rho_soil_avg)
    143 
    144 
    145 
    146 
    147 
    148 
     119  deallocate(rho_soil_avg)
    149120
    150121     do ig = 1,ngrid
    … …  
    164135        enddo
    165136      enddo
    166 
    167137   
    168138    deallocate(diff_rho)
    169 
    170 
    171 
    172 
    173 
    174139
    175140!=======================================================================
    … …  
    177142
    178143
    179 
    180 
    181 
    182144      END

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

    r2794 r2835  
    2323  REAL,    intent(in) ::  latitude(ngrid)          ! physical point field : Latitude
    2424  REAL,    intent(in) ::  initial_co2_ice(n_band_lat)  ! Initial/Actual surface of water ice
    25 
    26 
    2725
    2826!   OUTPUT
    … …  
    5957       present_co2(j).GT.initial_co2_ice(j)*(1+alpha_criterion)) then
    6058         STOPPING=.TRUE.
    61          print *, "j", j
    62          print *, "present_co2(j)", present_co2(j)
    63          print *, "initial_co2_ice(j)", initial_co2_ice(j)
    6459    endif
    6560  enddo

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

    r2794 r2835  
    2020!   INPUT
    2121  INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points
    22   REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
     22  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
    2323  REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
    2424  REAL,    intent(in) :: ini_surf
    … …  
    2727  REAL,    intent(in) :: global_ave_press_new
    2828
    29 
    3029!   OUTPUT
    3130  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
    … …  
    3332!   local:
    3433!   -----
    35   INTEGER :: i,islope                    ! Loop
     34  INTEGER :: i,islope   ! Loop
    3635  REAL :: present_surf  ! Initial/Actual surface of water ice
    3736
    … …  
    4645   do islope=1,nslope
    4746      if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
    48          print *, "i", i
    49          print *, "initial_h2o_ice(i,islope)", initial_h2o_ice(i,islope)
    50          print *, "qsurf(i,:)", qsurf(i,:)
    51          print *, "cell_area(i)", cell_area(i)
    52          print *, "present_surf",present_surf
    5347         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
    5448      endif
    5549   enddo
    5650  enddo
    57 
    58 !  print *, "initial_h2o_ice", initial_h2o_ice
    59 !  print *, "qsurf", qsurf
    60 
    61   print *, "present_surf", present_surf
    62   print *, "ini_surf", ini_surf
    63   print *, "ini_surf*0.8", ini_surf*(1-alpha_criterion)
    6451 
    6552!   check of the criterion
    6653  if(present_surf.LT.ini_surf*(1-alpha_criterion) .OR. &
    6754     present_surf.GT.ini_surf*(1+alpha_criterion)) then
    68   STOPPING=.TRUE.
     55    STOPPING=.TRUE.
     56    print *, "Reason of stopping : The surface of co2 ice sublimating reach the threshold:"
     57    print *, "Current surface of co2 ice sublimating=", present_surf
     58    print *, "Initial surface of co2 ice sublimating=", ini_surf
     59    print *, "Percentage of change accepted=", alpha_criterion*100
     60    print *, "present_surf<ini_surf*(1-alpha_criterion)", (present_surf.LT.ini_surf*(1-alpha_criterion))
    6961  endif
    7062
    … …  
    7365  endif
    7466
    75 !  if(global_ave_press_GCM.LT.global_ave_press_new*(1-alpha_criterion) .OR. &
    76 !     global_ave_press_GCM.GT.global_ave_press_new*(1+alpha_criterion)) then
    77 !  STOPPING=.TRUE.
    78 !  endif
    79 
    8067  if(global_ave_press_new.LT.global_ave_press_GCM*(0.9) .OR. &
    8168     global_ave_press_new.GT.global_ave_press_GCM*(1.1)) then
    82   STOPPING=.TRUE.
     69    STOPPING=.TRUE.
     70    print *, "Reason of stopping : The global pressure reach the threshold:"
     71    print *, "Current global pressure=", global_ave_press_new
     72    print *, "GCM global pressure=", global_ave_press_GCM
     73    print *, "Percentage of change accepted=", 0.1*100
     74    print *, "global_ave_press_new<global_ave_press_GCM*(0.9)", (global_ave_press_new.LT.global_ave_press_GCM*(0.9))
    8375  endif
    8476
    … …  
    8678
    8779
    88 
    89 
    90 

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

    r2794 r2835  
    2525  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
    2626
    27 
    2827!   OUTPUT
    2928  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
    … …  
    3938    STOPPING=.FALSE.
    4039
    41 !   computation of the actual surface
     40!   computation of the present surface of water ice sublimating
    4241  present_surf=0.
    4342  do i=1,ngrid
    4443    do islope=1, nslope
    4544      if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
    46 !         print *, "i", i
    47 !         print *, "initial_h2o_ice(i,islope)", initial_h2o_ice(i,islope)
    48 !         print *, "qsurf(i,islope)", qsurf(i,islope)
    49 !         print *, "cell_area(i)", cell_area(i)
    50 !         print *, "present_surf",present_surf
    5145         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
    5246      endif
    5347    enddo
    5448  enddo
    55 
    56 !  print *, "initial_h2o_ice", initial_h2o_ice
    57 !  print *, "qsurf", qsurf
    58 
    59 !  print *, "present_surf", present_surf
    60 !  print *, "ini_surf", ini_surf
    61 !  print *, "ini_surf*0.8", ini_surf*(1-alpha_criterion)
    6249 
    6350!   check of the criterion
    6451  if(present_surf.LT.ini_surf*(1-alpha_criterion) .OR. &
    6552     present_surf.GT.ini_surf*(1+alpha_criterion)) then
    66   STOPPING=.TRUE.
     53    STOPPING=.TRUE.
     54    print *, "Reason of stopping : The surface of water ice sublimating reach the threshold:"
     55    print *, "Current surface of water ice sublimating=", present_surf
     56    print *, "Initial surface of water ice sublimating=", ini_surf
     57    print *, "Percentage of change accepted=", alpha_criterion*100
     58    print *, "present_surf<ini_surf*(1-alpha_criterion)", (present_surf.LT.ini_surf*(1-alpha_criterion))
    6759  endif
    6860
    6961  if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
    70        STOPPING=.FALSE.
     62    STOPPING=.FALSE.
    7163  endif
    7264

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

    r2794 r2835  
    2121
    2222  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope   ! # of grid points along longitude/latitude/ total
    23 !  REAL, intent(in) ::  tendencies_h2o_ice_phys(ngrid) ! physical point field : Evolution of perenial ice over one year
    2423  REAL, intent(in) ::  cell_area(ngrid)
    2524
    … …  
    2827  LOGICAL :: STOPPING
    2928  REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    30 
    3129
    3230!   local:

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

    r2794 r2835  
    66
    77  USE temps_mod_evol, ONLY: dt_pem
    8           use comslope_mod, ONLY: subslope_dist
     8  use comslope_mod, ONLY: subslope_dist
    99
    1010      IMPLICIT NONE
    … …  
    2222
    2323  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope   ! # of grid points along longitude/latitude/ total
    24 !  REAL, intent(in) ::  tendencies_h2o_ice_phys(ngrid) ! physical point field : Evolution of perenial ice over one year
    2524  REAL, intent(in) ::  cell_area(ngrid)
    2625
    … …  
    3029  REAL, intent(inout) ::  tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    3130
    32 
    3331!   local:
    3432!   ----
    3533
    3634  INTEGER :: i,j,ig0,islope                                  ! loop variable
    37 !  REAL :: not_budget, budget
    3835  REAL :: pos_tend, neg_tend, real_coefficient,negative_part
    3936  REAL ::  new_tendencies(ngrid,nslope)
    4037
    41 
    4238!=======================================================================
    43 
    44 !  budget=sum(qsurf(:))
    4539
    4640  pos_tend=0.
    … …  
    5953  enddo
    6054
    61 !  print *, "pos_tend", pos_tend
    62 !  print *, "neg_tend", neg_tend
    63 
    6455  if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then
    6556     do i=1,ngrid
    6657       do islope=1,nslope
    6758       if(tendencies_h2o_ice_phys(i,islope).LT.0) then
    68 !          print *, "pos_tend/neg_tend", pos_tend/neg_tend
    6959          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(pos_tend/neg_tend)
    7060       else
    … …  
    7464     enddo
    7565  elseif(neg_tend.LT.pos_tend .and. neg_tend.GT.0) then
    76 !          print *, "neg_tend/pos_tend", neg_tend/pos_tend
    7766     do i=1,ngrid
    7867       do islope=1,nslope
    … …  
    8574     enddo
    8675  elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then
    87 
    88 !      call criterion_ice_stop(cell_area,1,qsurf*0.,STOPPING,ngrid,cell_area)
     76    print *, "Reason of stopping : There is either no water ice sublimating or no water ice increasing !!"
     77    print *, "Tendencies on ice sublimating=", neg_tend
     78    print *, "Tendencies on ice increasing=", pos_tend
     79    print *, "This can be due to the absence of water ice in the PCM run!!"
    8980      call criterion_ice_stop_water_slope(cell_area,1,qsurf(:,:)*0.,STOPPING,ngrid,cell_area)
    9081      do i=1,ngrid
    … …  
    9586  endif
    9687
    97  
    98  
    9988  negative_part = 0.
    100 
    10189
    10290! Evolution of the water ice for each physical point
    10391  do i=1,ngrid
    10492    do islope=1, nslope
    105 !      qsurf(i)=qsurf(i)+tendencies_h2o_ice_phys(i)*dt_pem
    10693      qsurf(i,islope)=qsurf(i,islope)+new_tendencies(i,islope)*dt_pem
    107 !      budget=budget+tendencies_h2o_ice_phys(i)*dt_pem
    10894      if (qsurf(i,islope).lt.0) then
    109 !        not_budget=not_budget+qsurf(i)
    110 !       print *, "NNqsurf(i,islope)", qsurf(i,islope)
    111 !       print *, "NNnew_tendencies(i,islope)", new_tendencies(i,islope)
    112 !       print *, "NNtendencies_h2o_ice_phys(i,islope)", tendencies_h2o_ice_phys(i,islope)
    11395        negative_part=negative_part-qsurf(i,islope)*cell_area(i)*subslope_dist(i,islope)
    11496        qsurf(i,islope)=0.
    11597        tendencies_h2o_ice_phys(i,islope)=0.
    116 !        print *, "NNineg", i
    117       endif
    118       if(qsurf(i,islope).NE.qsurf(i,islope)) then
    119 !          print *, "qsurf(i,islope)",qsurf(i,islope)
    120 !          print *, "new_tendencies",new_tendencies(i,islope)
    121 !          print *, "tendencies_h2o_ice_phys",tendencies_h2o_ice_phys(i,islope)
    122 !          print *, "i", i
    123 !          print *,"islope",islope
    12498      endif
    12599    enddo
    126100  enddo
    127101
    128 !  print *, "negative_part", negative_part
    129102  real_coefficient=negative_part/pos_tend
    130 !  print *, "real_coefficient", real_coefficient
    131103
    132104  do i=1,ngrid
    … …  
    140112
    141113
    142 ! Conservation of water ice
    143 !  qsurf(:)=qsurf(:)*budget/sum(qsurf(:))
    144 
    145 
    146114END SUBROUTINE evol_h2o_ice_s_slope

  • 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     

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

    r2794 r2835  
    5656  CALL err(NF90_CLOSE(fID),"close",fichnom)
    5757
     58  print *, "The number of timestep of the PCM run data=", timelen
     59
    5860  CONTAINS
    5961

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

    r2812 r2835  
    33
    44! I   Initialisation
    5 !    I_a READ run.def , run_pem.def
    6 !    I_b READ starfi_0.nc
    7 !    I_c READ GCM data and convert to the physical grid
    8 !    I_d Compute tendencies & Save initial situation
     5!    I_a READ run.def
     6!    I_b READ of start_evol.nc and starfi_evol.nc
     7!    I_c Subslope parametrisation
     8!    I_d READ GCM data and convert to the physical grid
     9!    I_e Initialisation of the PEM variable and soil
     10!    I_f Compute tendencies & Save initial situation
     11!    I_g Save initial PCM situation
     12!    I_h Read the PEMstart
     13!    I_i Compute orbit criterion
    914
    1015! II  Run
    1116!    II_a update pressure, ice and tracers
    12 !    II_b CO2 glaciers flows
    13 !    II_c Update surface and soil temperatures
    14 !    II_d Update the tendencies and stopping criterion
     17!    II_b Evolution of the ice
     18!    II_c CO2 glaciers flows
     19!    II_d Update surface and soil temperatures
     20!    II_e Update the tendencies
     21!    II_f Checking the stopping criterion
    1522
    1623! III Output
    17 !    III_a Recomp tendencies for the start and startfi
     24!    III_a Update surface value for the PCM start files
    1825!    III_b Write start and starfi.nc
    1926!    III_c Write start_pem
    2027
    2128!------------------------
    22 
    2329
    2430PROGRAM pem
    … …  
    3238                           hmons, summit, base,albedo_h2o_frost, &
    3339                          frost_albedo_threshold,emissiv
    34       use dimradmars_mod, only: totcloudfrac, albedo, &
    35                                 ini_dimradmars_mod
    36       use turb_mod, only: q2, wstar, ini_turb_mod
    37       use dust_param_mod, only: tauscaling, ini_dust_param_mod
    38 !      use co2cloud_mod, only: mem_Mccn_co2, mem_Mh2o_co2,&
    39 !      &                        mem_Nccn_co2,ini_co2cloud
     40      use dimradmars_mod, only: totcloudfrac, albedo
     41      use turb_mod, only: q2, wstar
     42      use dust_param_mod, only: tauscaling
    4043      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
    4144                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
    4245                        nf90_inquire_dimension,nf90_close
    4346      use phyredem, only: physdem0, physdem1
    44       use tracer_mod, only: noms,nqmx,igcm_h2o_ice ! tracer names
     47      use tracer_mod, only: noms,igcm_h2o_ice ! tracer names
    4548
    4649! For phyredem :
    … …  
    5659      USE iniphysiq_mod, ONLY: iniphysiq
    5760      USE infotrac
    58       USE temps_mod_evol, ONLY: nyear, dt_pem
    59 !     USE vertical_layers_mod, ONLY: ap,bp
    6061      USE comcstfi_h, only: pi
    6162
    … …  
    6869                           iflat
    6970
    70       USE geometry_mod, only: longitude_deg,latitude_deg
     71      USE geometry_mod, only: latitude_deg
    7172
    7273      use pemredem, only:  pemdem1
    73       USE soil_evolution_mod, ONLY: soil_pem_CN
     74
    7475!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
    75 
     76      USE soil_evolution_mod, ONLY: soil_pem, soil_pem_CN
    7677      use comsoil_h_PEM, only: ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
    77                               TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM, &        ! soil thermal inertia         
    78                               tsoil_PEM ,       &        !number of subsurface layers
    79                               mlayer_PEM,layer_PEM, &       ! soil mid layer depths
    80                               fluxgeo, &  ! geothermal flux
    81                               co2_adsorbded_phys ! mass of co2 in the regolith
     78                              TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM,         & ! soil thermal inertia         
     79                              tsoil_PEM, mlayer_PEM,layer_PEM,                  & !number of subsurface layers, soil mid layer depths
     80                              fluxgeo, co2_adsorbded_phys                         ! geothermal flux, mass of co2 in the regolith
    8281      use adsorption_mod, only : regolith_co2adsorption
     82
     83!!! For orbit parameters
     84      USE temps_mod_evol, ONLY: dt_pem, evol_orbit_pem, Max_iter_pem
     85      use planete_h, only: aphelie, periheli, year_day, peri_day, &
     86                          obliquit
     87      use orbit_param_criterion_mod, only : orbit_param_criterion
     88      use recomp_orb_param_mod, only: recomp_orb_param
     89
     90
    8391  IMPLICIT NONE
    8492
    … …  
    104112
    105113! Variable for reading start.nc
    106       character (len = *), parameter :: FILE_NAME_start = "start_0.nc" !Name of the file used for initialsing the PEM
     114      character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM
    107115  !   variables dynamiques
    108116  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    … …  
    120128! Variable for reading starfi.nc
    121129
    122       character (len = *), parameter :: FILE_NAME = "startfi_0.nc" !Name of the file used for initialsing the PEM
     130      character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
    123131      integer :: ncid, varid,status                                !Variable for handling opening of files
    124132      integer :: phydimid, subdimid, nlayerdimid, nqdimid          !Variable ID for Netcdf files
    … …  
    152160      REAL , dimension(:,:), allocatable :: min_h2o_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
    153161      REAL , dimension(:,:), allocatable :: min_h2o_ice_s_2        ! LON x LAT field : minimum of water ice at each point for the second year
    154       REAL , dimension(:,:,:), allocatable :: h2o_ice_first_last_day     ! LON x LAT x 2 field : First and Last value of a GCM year simulation of water ice
    155162
    156163      REAL , dimension(:,:), allocatable :: min_co2_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
    157164      REAL , dimension(:,:), allocatable :: min_co2_ice_s_2        ! LON x LAT field : minimum of water ice at each point for the second year
    158       REAL , dimension(:,:,:), allocatable :: co2_ice_first_last_day     ! LON x LAT x 2 field : First and Last value of a GCM year simulation of water ice
    159 
    160       REAL, dimension(:),allocatable  :: local_old_press           ! physical point field : Local pressure of initial/previous time step
    161       REAL, dimension(:),allocatable  :: local_new_press           ! physical point field : Local pressure of current time step
    162165
    163166      REAL :: global_ave_press_GCM
    … …  
    166169
    167170      REAL , dimension(:,:), allocatable ::  zplev_new
    168       REAL , dimension(:,:), allocatable :: zplev_old
     171      REAL , dimension(:,:), allocatable ::  zplev_gcm
    169172      REAL , dimension(:,:,:), allocatable ::  zplev_new_timeseries  ! same but with the time series
    170       REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries! same but with the time series
    171 
    172 
    173 
    174 
    175       REAL :: tot_co2_atm,tot_var_co2_atm
    176 
    177 
    178       INTEGER :: year_iter  !Counter for the number of PEM iteration
     173      REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries   ! same but with the time series
     174
    179175      LOGICAL :: STOPPING_water   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
    180176      LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
    181       LOGICAL :: STOPPING_co2   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
    182       LOGICAL :: STOPPING_1_co2 ! Logical : is there still water ice to sublimate?
     177      LOGICAL :: STOPPING_co2     ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
     178      LOGICAL :: STOPPING_1_co2   ! Logical : is there still water ice to sublimate?
    183179
    184180      REAL, dimension(:,:,:),allocatable  :: q_co2_GCM ! Initial amount of co2 in the first layer
    … …  
    195191      REAL, ALLOCATABLE ::  q_h2o_GCM(:,:,:)
    196192      REAL ,allocatable ::  q_h2o_PEM_phys(:,:)
    197       real ,allocatable :: q_phys(:,:,:)
    198193      integer :: timelen
    199194      REAL :: ave
    200195
    201196      REAL, ALLOCATABLE :: p(:,:)  !(ngrid,llmp1)
    202       REAL :: extra_mass ! Extra mass of a tracer if it is greater than 1
    203 
    204       REAL :: deficit_mass ! Extra mass of a tracer if it is lower than 0
    205 
    206 
    207 
     197      REAL :: extra_mass           ! Extra mass of a tracer if it is greater than 1
    208198
    209199!!!!!!!!!!!!!!!!!!!!!!!! SLOPE
    210200      character*2 :: str2
    211       REAL ,allocatable :: watercap_slope(:,:)    !(ngrid,nslope)
     201      REAL ,allocatable :: watercap_slope(:,:)                           !(ngrid,nslope)
    212202      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of water ice at each point for the first year
    213203      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of water ice at each point for the second year
    … …  
    217207
    218208      REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope        ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
    219       REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope        ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
    220 
    221 
    222       REAL, dimension(:,:),allocatable  :: initial_co2_ice_sublim_slope           ! physical point field : Logical array indicating sublimating point of co2 ice
     209      REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope     ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
     210
     211
     212      REAL, dimension(:,:),allocatable  :: initial_co2_ice_sublim_slope    ! physical point field : Logical array indicating sublimating point of co2 ice
    223213      REAL, dimension(:,:),allocatable  :: initial_h2o_ice_slope           ! physical point field : Logical array indicating sublimating point of h2o ice
    224214      REAL, dimension(:,:),allocatable  :: initial_co2_ice_slope           ! physical point field : Logical array indicating sublimating point of co2 ice
    … …  
    230220      REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical point field : Tendency of evolution of perenial co2 ice over a year
    231221
    232       REAL,SAVE,ALLOCATABLE,DIMENSION(:) ::  co2_hmax                          ! Maximum height  for CO2 deposit on slopes (m)
     222      REAL,SAVE,ALLOCATABLE,DIMENSION(:) ::  co2_hmax                      ! Maximum height  for CO2 deposit on slopes (m)
    233223
    234224      REAL, PARAMETER :: rho_co2 = 1600           ! CO2 ice density (kg/m^3)
    … …  
    237227      REAL , dimension(:), allocatable :: flag_co2flow_mesh(:)  !(ngrid)          ! To flag where there is a glacier flow
    238228
    239 
    240229!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL
    241230
    242 
    243      REAL, ALLOCATABLE :: tsurf_ave(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
    244      REAL, ALLOCATABLE  :: tsurf_ave_phys(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature [K]
    245      REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
    246      REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
     231     REAL, ALLOCATABLE :: tsurf_ave(:,:,:)          ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
     232     REAL, ALLOCATABLE  :: tsurf_ave_phys(:,:)      ! IG x LAT x SLOPE field : Averaged Surface Temperature [K]
     233     REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:)        ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
     234     REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:)    ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
    247235     REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) !IG x SLOPE field : Averaged Soil Temperature [K]
    248236
    249      REAL, ALLOCATABLE :: TI_GCM(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Thermal Inertia  [SI]
    250      REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:) !LONX LAT x SLOPE XTULES field : NOn averaged Surf Temperature [K]
     237     REAL, ALLOCATABLE :: TI_GCM(:,:,:,:)                  ! LON x LAT x SLOPE field : Averaged Thermal Inertia  [SI]
     238     REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:)    !LONX LAT x SLOPE XTULES field : NOn averaged Surf Temperature [K]
    251239     REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) !IG x SLOPE XTULES field : NOn averaged Surf Temperature [K]
    252240
    253241     REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    254      REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    255 
    256      REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
     242     REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:)    !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
     243
     244     REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:)     ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
    257245     REAL, ALLOCATABLE  :: tsurf_ave_phys_yr1(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
    258246
    259 
    260 
    261      REAL :: tol_soil = 1.e-6
    262      INTEGER :: countmax = 15
    263      INTEGER :: countloop
    264      REAL, ALLOCATABLE :: tsoil_phys_PEM_saved(:,:,:)
    265      REAL :: error_tsoil
    266 
    267 
    268      LOGICAL :: firstcall
    269 !    REAL, PARAMETER :: daysec=88775.              ! duree du sol (s)  ~88775 s
    270      REAL, PARAMETER :: year_day = 669
    271247     REAL, PARAMETER :: year_step = 1
     248     INTEGER :: year_iter        !Counter for the number of PEM iteration
     249     INTEGER :: year_iter_max
    272250     REAL :: timestep
    273251
    274      REAL, ALLOCATABLE :: inertiesoil(:,:) ! Thermal inertia of the mesh for restart
    275 
    276      REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:) ! Averaged GCM Thermal Inertia  [SI]
     252     REAL, ALLOCATABLE :: inertiesoil(:,:)    ! Thermal inertia of the mesh for restart
     253
     254     REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:)  ! Averaged GCM Thermal Inertia  [SI]
    277255     REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Averaged GCM Thermal Inertia  [SI]
    278 
    279      REAL,ALLOCATABLE  :: q_h2o_PEM_phys_ave(:) ! averaged water vapor content
    280      REAL,ALLOCATABLE  :: interp_coef(:)
    281256
    282257     REAL,ALLOCATABLE  :: ice_depth(:,:)
    283258     REAL,ALLOCATABLE  :: TI_locslope(:,:)
    284259     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:)
    285      REAL,ALLOCATABLE  :: Tsoil_locslope_CN(:,:)
    286260     REAL,ALLOCATABLE  :: Tsurf_locslope(:)
    287      REAL,ALLOCATABLE  :: Tsurf_locslope_prev(:)
    288261     REAL,ALLOCATABLE  :: alph_locslope(:,:)
    289262     REAL,ALLOCATABLE  :: beta_locslope(:,:)   
    290      REAL :: kcond
     263
    291264     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)
    292265     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)
    293266     REAL :: totmass_adsorbded
     267
    294268!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    295269
    … …  
    310284      time_phys=0. !test
    311285
    312 !      n_band_lat=18
     286! Some constants
     287
     288      ngrid=ngridmx
     289      nlayer=llm
    313290
    314291      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
    … …  
    317294      B=1/m_noco2
    318295
    319 
     296      year_day=669
     297      daysec=88775.
     298      dtphys=0
     299      timestep=year_day*daysec/year_step
    320300
    321301!------------------------
    … …  
    324304!    I_a READ run.def
    325305
    326 
    327306!------------------------
    328307
    329308!----------------------------READ run.def ---------------------
    330309      CALL conf_gcm( 99, .TRUE. )
    331       CALL conf_evol( 99, .TRUE. )
    332 
     310
     311      call infotrac_init
     312      nq=nqtot
    333313
    334314!------------------------
    … …  
    336316! I   Initialisation
    337317!    I_a READ run.def
    338 !    I_b READ starfi_0.nc
    339 
    340 !----------------------------READ startfi.nc ---------------------
    341 
    342 !----------------------------initialisation ---------------------
    343 
     318!    I_b READ of start_evol.nc and starfi_evol.nc
     319
     320!------------------------
     321
     322!----------------------------Initialisation : READ some constant of startfi_evol.nc ---------------------
     323
     324! In the gcm, these values are given to the physic by the dynamic.
     325! Here we simply read them in the startfi_evol.nc file
    344326      status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
    345       status =nf90_inq_dimid(ncid, "physical_points", phydimid)
    346       status =nf90_inquire_dimension(ncid, phydimid, len = ngrid)
    347 
    348       status =nf90_inq_dimid(ncid, "nlayer", nlayerdimid)
    349       status =nf90_inquire_dimension(ncid, nlayerdimid, len = nlayer)
    350 
    351       status =nf90_inq_dimid(ncid,"number_of_advected_fields",nqdimid)
    352       status =nf90_inquire_dimension(ncid, nqdimid, len = nq)
    353327
    354328      allocate(longitude(ngrid))
    … …  
    372346      status =nf90_close(ncid)
    373347
    374 
    375 
    376348!----------------------------READ start.nc ---------------------
    377    
    378      call infotrac_init
    379349
    380350     allocate(q(ip1jmp1,llm,nqtot))
    … …  
    392362     status = nf90_get_var(ncid, bpvarid, bp)
    393363     status =nf90_close(ncid)
    394      
    395        
    396    
    397     daysec=88775.
    398     dtphys=0
    399 
    400  
    401 
    402     CALL iniphysiq(iim,jjm,llm, &
     364
     365     CALL iniphysiq(iim,jjm,llm, &
    403366          (jjm-1)*iim+2,comm_lmdz, &
    404367          daysec,day_ini,dtphys/nsplit_phys, &
    405368          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
    406369          iflag_phys)
    407    DO nnq=1,nqtot
    408     if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
    409    ENDDO
    410 
    411  
    412 !----------------------------reading ---------------------
    413 
     370
     371     DO nnq=1,nqtot
     372       if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
     373     ENDDO
     374
     375!----------------------------READ startfi.nc ---------------------
    414376
    415377! First we read the initial state (starfi.nc)
    416378
    417 
    418        allocate(watercap_slope(ngrid,nslope))
    419        allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
    420        allocate(q_h2o_PEM_phys_ave(ngrid))
     379      allocate(watercap_slope(ngrid,nslope))
     380      allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
    421381      allocate(inertiesoil(ngrid,nsoilmx))
    422382
    423 
    424          CALL phyetat0 (FILE_NAME,0,0, &
     383      CALL phyetat0 (FILE_NAME,0,0, &
    425384              nsoilmx,ngrid,nlayer,nq,   &
    426385              day_ini,time_phys,         &
    … …  
    432391              qsurf_slope,watercap_slope)
    433392
    434 
    435 
    436      deallocate(TI_GCM_start) !not used then
    437 
     393     if(soil_pem) then
     394       deallocate(TI_GCM_start) !not used then
     395     endif
     396
     397! Remove unphysical values of surface tracer
    438398     DO i=1,ngrid
    439          DO nnq=1,nqtot
    440            DO islope=1,nslope
    441              if(qsurf_slope(i,nnq,islope).LT.0) then
    442                 qsurf_slope(i,nnq,islope)=0.
    443              endif
    444            enddo
     399       DO nnq=1,nqtot
     400         DO islope=1,nslope
     401           if(qsurf_slope(i,nnq,islope).LT.0) then
     402             qsurf_slope(i,nnq,islope)=0.
     403           endif
    445404         enddo
    446405       enddo
     406     enddo
     407
     408!------------------------
     409
     410! I   Initialisation
     411!    I_a READ run.def
     412!    I_b READ of start_evol.nc and starfi_evol.nc
     413!    I_c Subslope parametrisation
     414
     415!------------------------
     416
     417!----------------------------Subslope parametrisation definition ---------------------
    447418
    448419!     Define some slope statistics
    449 
    450 
    451      
    452 
    453        iflat=1
    454        DO islope=2,nslope
    455          IF(abs(def_slope_mean(islope)).lt. &
    456            abs(def_slope_mean(iflat)))THEN
    457            iflat = islope
    458          ENDIF     
    459 
    460        ENDDO
    461        PRINT*,'Flat slope for islope = ',iflat
    462        PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
     420     iflat=1
     421     DO islope=2,nslope
     422       IF(abs(def_slope_mean(islope)).lt. &
     423         abs(def_slope_mean(iflat))) THEN
     424         iflat = islope
     425       ENDIF     
     426     ENDDO
     427
     428     PRINT*,'Flat slope for islope = ',iflat
     429     PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
     430
    463431! CO2 max thickness (for glaciers flows)
    464 
    465        allocate(co2_hmax(nslope)
    466        if(nslope.eq.7) ! ugly way to implement that ...
     432       allocate(co2_hmax(nslope))
     433       if(nslope.eq.7) then ! ugly way to implement that ...
    467434!        CF documentation that explain how the values are computed
    468435         co2_hmax(1) = 1.5
    469          co2_hmax(7) = co2_max(1)
     436         co2_hmax(7) = co2_hmax(1)
    470437         co2_hmax(2) = 2.4
    471438         co2_hmax(6) = co2_hmax(2)
    … …  
    475442       endif
    476443
    477 
    478 
    479 
     444     allocate(flag_co2flow(ngrid,nslope))
     445     allocate(flag_co2flow_mesh(ngrid))
    480446
    481447       flag_co2flow(:,:) = 0.     
    482448       flag_co2flow_mesh(:) = 0.
    483449
    484 !------------------------
     450!---------------------------- READ GCM data ---------------------
    485451
    486452! I   Initialisation
    487 !    I_a READ run.def , run_pem.def
    488 !    I_b READ starfi_0.nc
    489 !    I_c READ GCM data and convert to the physical grid
     453!    I_a READ run.def
     454!    I_b READ of start_evol.nc and starfi_evol.nc
     455!    I_c Subslope parametrisation
     456!    I_d READ GCM data and convert to the physical grid
    490457
    491458!------------------------
    … …  
    493460! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the GCM run, saving only the minimum value
    494461
    495 
    496      call nb_time_step_GCM("concat_year_one.nc",timelen)
     462     call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
     463
    497464     allocate(min_h2o_ice_s_1(iim+1,jjm+1))
    498465     allocate(min_co2_ice_s_1(iim+1,jjm+1))
    … …  
    506473     allocate(tsurf_ave(iim+1,jjm+1,nslope))
    507474     allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope))
    508      allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
    509      allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
    510      allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
    511475     allocate(tsurf_ave_phys(ngrid,nslope))
    512476     allocate(tsurf_ave_phys_yr1(ngrid,nslope))
    513      allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
    514      allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
    515      allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    516 
    517477     allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen))
    518478     allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen))
    … …  
    520480     allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen))
    521481     allocate(Tsurfave_before_saved(ngrid,nslope))
    522      allocate(tsoil_phys_PEM_saved(ngrid,nsoilmx_PEM,nslope))
     482     allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
     483     allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
     484     allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
     485     allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
     486     allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
     487     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    523488     allocate(delta_co2_adsorbded(ngrid))
    524      allocate(q_phys(ngrid,llm,nqtot))
    525 
    526 
    527 
    528      call read_data_GCM("concat_year_one.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,&   
     489
     490     print *, "Downloading data Y1..."
     491
     492     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,&   
    529493                       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)
    530494
    531 
    532495! 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
     496
     497     print *, "Downloading data Y1 done"
    533498
    534499     allocate(min_h2o_ice_s_2(iim+1,jjm+1))
    … …  
    537502     allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
    538503
    539      call read_data_GCM("concat_year_two.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, &
     504     print *, "Downloading data Y2"
     505
     506     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, &
    540507                  nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope)
    541508
    542 
    543 
     509     print *, "Downloading data Y2 done"
    544510
    545511! The variables in the dynamic grid are transfered to the physical grid
    546 
    547512
    548513     allocate(vmr_co2_gcm_phys(ngrid,timelen))
    … …  
    552517     allocate(q_co2_GCM_phys(ngrid,timelen))
    553518     allocate(q_co2_PEM_phys(ngrid,timelen))
    554 
     519     allocate(ps_phys(ngrid))
     520     allocate(ps_phys_timeseries(ngrid,timelen))
     521     allocate(ps_phys_timeseries_yr1(ngrid,timelen))
    555522
    556523     CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys)
    557524     CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys)
    558525     CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys)
    559 
    560      do nnq=1,nqtot
    561 
    562         call gr_dyn_fi(llm,iip1,jjp1,ngridmx,q(:,:,nnq),q_phys(:,:,nnq))
    563      enddo
    564 
     526     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
     527     call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
     528     call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
     529     CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
     530     CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
    565531
    566532     deallocate(vmr_co2_gcm)
    567533     deallocate(q_h2o_GCM)
    568534     deallocate(q_co2_GCM)
    569 
    570 
    571     q_co2_PEM_phys(:,:)=  q_co2_GCM_phys(:,:)
    572     q_h2o_PEM_phys(:,:)=  q_h2o_GCM_phys(:,:)
    573 
    574 
    575      allocate(ps_phys(ngrid))
    576      allocate(ps_phys_timeseries(ngrid,timelen))
    577      allocate(ps_phys_timeseries_yr1(ngrid,timelen))
    578 
    579 
    580 
    581 
    582 
    583      call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
    584 
    585 
    586 ! for the time series:
    587 
    588 
    589 
    590     call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
    591     call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
    592     deallocate(ps_GCM)
    593     deallocate(ps_GCM_yr1)
    594 
    595 
    596 !!!!!!!!!!!!!!!!!!!!! Initialisation for soil_PEM
     535     deallocate(ps_GCM)
     536     deallocate(ps_GCM_yr1)
     537     deallocate(tsurf_ave)
     538     deallocate(tsurf_ave_yr1)
     539
     540     q_co2_PEM_phys(:,:)=  q_co2_GCM_phys(:,:)
     541     q_h2o_PEM_phys(:,:)=  q_h2o_GCM_phys(:,:)
     542
     543!------------------------
     544
     545! I   Initialisation
     546!    I_a READ run.def
     547!    I_b READ of start_evol.nc and starfi_evol.nc
     548!    I_c Subslope parametrisation
     549!    I_d READ GCM data and convert to the physical grid
     550!    I_e Initialisation of the PEM variable and soil
     551
     552!------------------------
     553
     554!---------------------------- Initialisation of the PEM soil and values ---------------------
    597555
    598556      call end_comsoil_h_PEM
    599557      call ini_comsoil_h_PEM(ngrid,nslope)
    600       timestep=year_day*daysec/year_step
    601 
    602558
    603559      allocate(ice_depth(ngrid,nslope))
    … …  
    605561
    606562      DO islope = 1,nslope
     563        if(soil_pem) then
    607564          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope))
     565        endif !soil_pem
     566        DO t=1,timelen
     567          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
     568          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
     569        enddo
    608570      ENDDO
    609571
    610 
     572      deallocate(co2_ice_GCM_slope)
     573      deallocate(TI_GCM)
     574
     575  if(soil_pem) then
     576
     577      deallocate(tsurf_GCM_timeseries)
    611578      call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
    612579
    613       deallocate(TI_GCM)
    614 
    615 
    616580      DO islope = 1,nslope
    617 
    618581        DO l=1,nsoilmx
    619582           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
     583           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
     584           DO t=1,timelen
     585             CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
     586           ENDDO
    620587        ENDDO
    621588        DO l=nsoilmx+1,nsoilmx_PEM
    622589           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
    623         ENDDO
    624       ENDDO
    625 
    626 
    627 
    628        deallocate(tsoil_ave_yr1)
    629 
    630 
    631 
    632 
    633 
    634 
    635 
    636 
    637 
    638 
    639       DO islope = 1,nslope
    640 
    641         DO l=1,nsoilmx
    642            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
    643         ENDDO
    644         DO l=nsoilmx+1,nsoilmx_PEM
    645590           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
    646591        ENDDO
    647592      ENDDO
    648      
    649 
    650 
    651        deallocate(tsoil_ave)
    652 
    653 
    654       DO islope = 1,nslope
    655        
    656         DO l=1,nsoilmx
    657           DO t=1,timelen
    658            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
    659           ENDDO
    660         ENDDO
    661       ENDDO
    662 
     593
     594      deallocate(tsoil_ave_yr1)
     595      deallocate(tsoil_ave)
    663596      deallocate(tsoil_GCM_timeseries)
    664      
    665  
    666 
    667 
    668        CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
    669        
    670        CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
    671 
    672        deallocate(tsurf_ave)
    673        deallocate(tsurf_ave_yr1)
    674 
    675    
    676 
    677 
    678       DO islope = 1,nslope
    679 
    680         DO t=1,timelen
    681           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
    682         enddo
    683       enddo
    684       deallocate(tsurf_GCM_timeseries)
    685 
    686 
    687   DO islope = 1,nslope
    688 
    689         DO t=1,timelen
    690           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
    691         enddo
    692       enddo
    693 
    694 
    695 
    696       deallocate(co2_ice_GCM_slope)
    697 
    698 
     597  endif !soil_pem
    699598
    700599!------------------------
    701600
    702601! I   Initialisation
    703 !    I_a READ run.def , run_pem.def
    704 !    I_b READ starfi_0.nc
    705 !    I_c READ GCM data and convert to the physical grid
    706 !    I_d Compute tendencies & Save initial situation
    707 
    708 !------------------------
    709 
    710 
    711 
    712 
    713 
    714 !----- Compute tendencies
    715 
     602!    I_a READ run.def
     603!    I_b READ of start_evol.nc and starfi_evol.nc
     604!    I_c Subslope parametrisation
     605!    I_d READ GCM data and convert to the physical grid
     606!    I_e Initialisation of the PEM variable and soil
     607!    I_f Compute tendencies & Save initial situation
     608
     609!----- Compute tendencies from the PCM run
    716610
    717611     allocate(tendencies_h2o_ice(iim+1,jjm+1))
    718612     allocate(tendencies_h2o_ice_phys(ngrid))
    719 
    720613     allocate(tendencies_co2_ice(iim+1,jjm+1))
    721614     allocate(tendencies_co2_ice_phys(ngrid))
    722 
    723615     allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
    724616     allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
    725617     allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope))
    726 
    727 
    728618     allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope))
    729619     allocate(tendencies_h2o_ice_phys_slope(ngrid,nslope))
    730620
    731 !  Compute the tendencies of the evolution of water ice over the years
     621!  Compute the tendencies of the evolution of ice over the years
    732622
    733623      call compute_tendencies(tendencies_h2o_ice,min_h2o_ice_s_1,&
    … …  
    740630             min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope)
    741631
    742 
    743632      tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:)
    744633
    … …  
    746635             min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope)
    747636
    748      
    749 !----- Save initial situation
     637!------------------------
     638
     639! I   Initialisation
     640!    I_a READ run.def
     641!    I_b READ of start_evol.nc and starfi_evol.nc
     642!    I_c Subslope parametrisation
     643!    I_d READ GCM data and convert to the physical grid
     644!    I_e Initialisation of the PEM variable and soil
     645!    I_f Compute tendencies & Save initial situation
     646!    I_g Save initial PCM situation
     647
     648!---------------------------- Save initial PCM situation ---------------------
    750649
    751650     allocate(initial_h2o_ice(ngrid))
    … …  
    754653     allocate(initial_co2_ice_slope(ngrid,nslope))
    755654     allocate(initial_h2o_ice_slope(ngrid,nslope))
    756      year_iter=0
    757655
    758656! We save the places where water ice is sublimating
     657! We compute the surface of water ice sublimating
     658  ini_surf=0.
     659  ini_surf_co2=0.
     660  ini_surf_h2o=0.
     661  Total_surface=0.
    759662  do i=1,ngrid
    760       if (tendencies_h2o_ice_phys(i).LT.0) then
     663    Total_surface=Total_surface+cell_area(i)
     664    if (tendencies_h2o_ice_phys(i).LT.0) then
    761665         initial_h2o_ice(i)=1.
    762       else
     666         ini_surf=ini_surf+cell_area(i)
     667    else
    763668         initial_h2o_ice(i)=0.         
    764       endif
     669    endif
    765670    do islope=1,nslope
    766 
    767671      if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
    768672         initial_co2_ice_sublim_slope(i,islope)=1.
     673         ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
    769674      else
    770675         initial_co2_ice_sublim_slope(i,islope)=0.         
    … …  
    777682      if (tendencies_h2o_ice_phys_slope(i,islope).LT.0) then
    778683         initial_h2o_ice_slope(i,islope)=1.
     684         ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
    779685      else
    780686         initial_h2o_ice_slope(i,islope)=0.         
    … …  
    783689  enddo
    784690
    785 ! We compute the surface of water ice sublimating
    786   ini_surf=0.
    787   ini_surf_co2=0.
    788   ini_surf_h2o=0.
    789   Total_surface=0.
    790   do i=1,ngrid
    791     if (initial_h2o_ice(i).GT.0.5) then
    792        ini_surf=ini_surf+cell_area(i)
    793     endif
    794     do islope=1,nslope
    795       if (initial_co2_ice_sublim_slope(i,islope).GT.0.5) then
    796          ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
    797       endif
    798       if (initial_h2o_ice_slope(i,islope).GT.0.5) then
    799          ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
    800       endif
    801     enddo
    802     Total_surface=Total_surface+cell_area(i)
    803   enddo
    804 
    805      print *, "ini_surf_co2=", ini_surf_co2
    806      print *, "ini_surf=", ini_surf
    807      print *, "ini_surf_h2o=", ini_surf_h2o
    808 
    809 
    810      allocate(local_old_press(ngrid))
    811      allocate(local_new_press(ngrid))
     691     print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
     692     print *, "Total initial surface of h2o ice sublimating=", ini_surf
     693     print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
     694     print *, "Total surface of the planet=", Total_surface
    812695   
    813      allocate(zplev_new(ngrid,nlayer+1))
    814      allocate(zplev_old(ngrid,nlayer+1))
    815 
     696     allocate(zplev_gcm(ngrid,nlayer+1))
     697
     698     DO l=1,nlayer+1
     699       DO ig=1,ngrid
     700         zplev_gcm(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
     701       ENDDO
     702     ENDDO
    816703
    817704     global_ave_press_old=0.
    … …  
    821708
    822709     global_ave_press_GCM=global_ave_press_old
    823      print *, "global_ave_press_old", global_ave_press_old
    824      
    825 
    826 
    827 !----- Time loop
    828      print *, "Before timeloop"
    829      allocate(flag_co2flow(ngrid,nslope))
    830      allocate(flag_co2flow_mesh(ngrid))
    831 
    832      firstcall=.TRUE.
    833 
    834 
    835 
    836 
     710     global_ave_press_new=global_ave_press_old
     711     print *, "Initial global average pressure=", global_ave_press_GCM
    837712
    838713!------------------------
    839714
    840715! I   Initialisation
    841 !    I_a READ run.def , run_pem.def
    842 !    I_b READ starfi_0.nc
    843 !    I_c READ GCM data and convert to the physical grid
    844 !    I_d Compute tendencies & Save initial situation
    845 !    I_e Read startfi_PEM
    846 
    847 !------------------------
    848 
    849 
    850 ! now we complete with the PEMstart
     716!    I_a READ run.def
     717!    I_b READ of start_evol.nc and starfi_evol.nc
     718!    I_c Subslope parametrisation
     719!    I_d READ GCM data and convert to the physical grid
     720!    I_e Initialisation of the PEM variable and soil
     721!    I_f Compute tendencies & Save initial situation
     722!    I_g Save initial PCM situation
     723!    I_h Read the PEMstart
     724
     725!---------------------------- Read the PEMstart ---------------------
    851726
    852727      call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
    853       tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries_yr1,ps_phys_timeseries,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,&
     728      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,&
    854729      tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), co2_adsorbded_phys,delta_co2_adsorbded)
    855 totmass_adsorbded = 0.
     730
     731
     732    if(soil_pem) then
     733     totmass_adsorbded = 0.
     734
    856735     do ig = 1,ngrid
    857736      do islope =1, nslope
    858737       do l = 1,nsoilmx_PEM - 1
    859 
    860738          totmass_adsorbded = totmass_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
    861739          subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
    862740          cell_area(ig)
    863741       enddo
    864 
    865742      enddo
    866 
    867743     enddo
    868744
    869      write(*,*) "tot mass in the regolith", totmass_adsorbded
    870      stop
     745     write(*,*) "Tot mass in the regolith=", totmass_adsorbded
     746     deallocate(tsoil_ave_phys_yr1) 
     747    endif !soil_pem
    871748     deallocate(tsurf_ave_phys_yr1)
    872      deallocate(ps_phys_timeseries_yr1)
    873      deallocate(tsoil_ave_phys_yr1)   
    874 
     749     deallocate(ps_phys_timeseries_yr1) 
     750
     751!------------------------
     752
     753! I   Initialisation
     754!    I_a READ run.def
     755!    I_b READ of start_evol.nc and starfi_evol.nc
     756!    I_c Subslope parametrisation
     757!    I_d READ GCM data and convert to the physical grid
     758!    I_e Initialisation of the PEM variable and soil
     759!    I_f Compute tendencies & Save initial situation
     760!    I_g Save initial PCM situation
     761!    I_h Read the PEMstar
     762!    I_i Compute orbit criterion
     763
     764     CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
     765
     766     if(evol_orbit_pem) then
     767       call orbit_param_criterion(year_iter_max)
     768     else
     769       year_iter_max=Max_iter_pem
     770     endif
    875771
    876772!--------------------------- END INITIALISATION ---------------------
    877773
    878 
    879 
    880 
    881774!---------------------------- RUN ---------------------
    882 
    883 
    884 
    885775
    886776!------------------------
    … …  
    890780
    891781!------------------------
    892 
    893    
    894      do while (.true.)
    895 ! II.a.1. global pressure
    896      global_ave_press_new=global_ave_press_old
    897 
     782     year_iter=0
     783     print *, "Max_iter_pem", Max_iter_pem
     784     do while (year_iter.LT.year_iter_max)
     785
     786! II.a.1. Compute updated global pressure
     787     print *, "Recomputing the new pressure..."
    898788     do i=1,ngrid
    899789       do islope=1,nslope
    900790           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
    901791       enddo
    902      
    903       global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
     792       if(soil_pem) then
     793         global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
     794       endif
    904795     enddo
    905 ! II.a.2. tracers
     796
     797! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
    906798     allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
    907         DO l=1,nlayer+1
    908          DO ig=1,ngrid
    909           zplev_old(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
    910           zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
    911          ENDDO
    912         ENDDO
    913 
    914      do i=1,ip1jmp1
    915 
    916        ps(i)=ps(i)*global_ave_press_new/global_ave_press_old
    917 
    918      enddo
     799     print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
     800     DO l=1,nlayer+1
     801       DO ig=1,ngrid
     802         zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
     803       ENDDO
     804     ENDDO
     805
     806! II.a.3. Surface pressure timeseries
     807     print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
    919808     do i = 1,ngrid
    920        ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_old
    921809       ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old
    922810     enddo
    923811
     812! II.a.4. New pressure levels timeseries
    924813     allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
    925 
    926         do l=1,nlayer+1
    927          do ig=1,ngrid
    928           zplev_new(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
    929           zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
    930          enddo
    931         enddo
    932 
    933 
    934       print *, "1 is okay"
    935 
    936         DO nnq=1,nqtot
    937         if (noms(nnq).NE."co2") then
    938           DO l=1,llm-1
    939             DO ig=1,ngrid
    940                q(ig,l,nnq)=q(ig,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))
    941             ENDDO
    942             q(:,llm,nnq)=q(:,llm-1,nnq)
    943           ENDDO
    944         else
    945           DO l=1,llm-1
    946             DO ig=1,ngrid
    947                 q(ig,l,nnq)=q(ig,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))  &
    948                                 + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
    949                                       (zplev_old(ig,l)-zplev_old(ig,l+1))     )  / &
    950                                       (zplev_new(ig,l)-zplev_new(ig,l+1))
    951             ENDDO
    952            q(:,llm,nnq)=q(:,llm-1,nnq)
    953           ENDDO
    954         endif
    955         ENDDO
    956 
    957 
    958  DO nnq=1,nqtot
    959        DO ig=1,ngrid
    960          DO l=1,llm-1
    961           if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number" .or. noms(nnq).NE."ccn_number") ) then
    962               extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
    963               q(ig,l,nnq)=1.
    964               q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
    965           endif
    966           if(q(ig,l,nnq).LT.0) then
    967               q(ig,l,nnq)=1E-30
    968           endif
    969          ENDDO
     814     print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
     815     do l=1,nlayer+1
     816       do ig=1,ngrid
     817         zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
    970818       enddo
    971    enddo
    972 
    973 
    974 
    975 
    976   l=1
    977   DO ig=1,ngrid
    978       DO t = 1, timelen
    979 
     819     enddo
     820
     821! II.a.5. Tracers timeseries
     822      print *, "Recomputing of tracer VMR timeseries for the new pressure..."
     823
     824     l=1
     825     DO ig=1,ngrid
     826       DO t = 1, timelen
    980827         q_h2o_PEM_phys(ig,t)=q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
    981828         if(q_h2o_PEM_phys(ig,t).LT.0) then
    … …  
    985832           q_h2o_PEM_phys(ig,t)=1.
    986833         endif
    987 
    988    enddo
    989   enddo
    990 
    991 
    992 
    993   DO ig=1,ngrid
    994       DO t = 1, timelen
     834       enddo
     835     enddo
     836
     837     DO ig=1,ngrid
     838       DO t = 1, timelen
    995839         q_co2_PEM_phys(ig,t)=q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))  &
    996840                                + (   (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))  -       &
    … …  
    1004848         mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
    1005849         vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
    1006       ENDDO
    1007   ENDDO
     850       ENDDO
     851     ENDDO
    1008852   
    1009853     deallocate(zplev_new_timeseries)
    1010854     deallocate(zplev_old_timeseries)
    1011855
    1012 
    1013 ! II.a.3. Evolution of the ice
    1014 
    1015 
    1016      call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
    1017       print *, "4 is okay"
    1018 
    1019      call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
    1020       print *, "5 is okay"
    1021 
    1022 
    1023 !------------------------
    1024 
    1025856! II  Run
    1026857!    II_a update pressure, ice and tracers
    1027 !    II_b CO2 glaciers flows
    1028 
    1029 !------------------------
    1030 
    1031    
    1032 
    1033 
     858!    II_b Evolution of the ice
     859
     860! II.b. Evolution of the ice
     861      print *, "Evolution of h2o ice"
     862     call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
     863
     864      print *, "Evolution of co2 ice"
     865     call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
     866
     867!------------------------
     868
     869! II  Run
     870!    II_a update pressure, ice and tracers
     871!    II_b Evolution of the ice
     872!    II_c CO2 glaciers flows
     873
     874!------------------------
     875
     876      print *, "Co2 glacier flow"
    1034877       DO ig = 1,ngrid
    1035878        DO islope = 1,nslope
    … …  
    1038881            IF(co2ice_slope(ig,islope).ge.rho_co2*co2_hmax(islope) * &
    1039882                  cos(pi*def_slope_mean(islope)/180.)) THEN
    1040 
    1041883! Second: determine the flatest slopes possible:
    1042        
    1043884                IF(islope.gt.iflat) THEN
    1044885                  iaval=islope-1
    … …  
    1048889                do while ((iaval.ne.iflat).and.  &
    1049890                    (subslope_dist(ig,iaval).eq.0))
    1050 
    1051891                  IF(iaval.gt.iflat) THEN
    1052892                     iaval=iaval-1
    … …  
    1054894                     iaval=iaval+1
    1055895                  ENDIF
    1056              
    1057896                enddo
    1058 
    1059 
    1060                 co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + &
     897              co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + &
    1061898               (co2ice_slope(ig,islope) - rho_co2* co2_hmax(islope) *     &
    1062899               cos(pi*def_slope_mean(islope)/180.)) *             &
    … …  
    1065902               cos(pi*def_slope_mean(islope)/180.)               
    1066903
    1067                 co2ice_slope(ig,islope)=rho_co2*co2_hmax(islope) *        &
    1068                   cos(pi*def_slope_mean(islope)/180.)
    1069 
    1070              flag_co2flow(ig,islope) = 1.
    1071              flag_co2flow_mesh(ig) = 1.
     904              co2ice_slope(ig,islope)=rho_co2*co2_hmax(islope) *        &
     905               cos(pi*def_slope_mean(islope)/180.)
     906
     907              flag_co2flow(ig,islope) = 1.
     908              flag_co2flow_mesh(ig) = 1.
    1072909            ENDIF ! co2ice > hmax
    1073910          ENDIF ! iflattsoil_lope
    … …  
    1075912       ENDDO !ig
    1076913
    1077       print *, "6 is okay"
    1078 
    1079914!------------------------
    1080915
    1081916! II  Run
    1082917!    II_a update pressure, ice and tracers
    1083 !    II_b CO2 glaciers flows
    1084 !    II_c Update surface and soil temperatures
    1085 
    1086 !------------------------
    1087 
    1088 
    1089 ! II_c.1 Update Tsurf
     918!    II_b Evolution of the ice
     919!    II_c CO2 glaciers flows
     920!    II_d Update surface and soil temperatures
     921
     922!------------------------
     923
     924! II_d.1 Update Tsurf
     925
     926      print *, "Updating the new Tsurf"
    1090927      bool_sublim=0
    1091928      Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:)
    1092        DO ig = 1,ngrid
     929      DO ig = 1,ngrid
    1093930        DO islope = 1,nslope
    1094            if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-5) THEN !co2ice disappeared, look for closest point without co2ice
    1095 
    1096 
    1097               if(latitude_deg(ig).gt.0) then
    1098                 do ig_loop=ig,ngrid
    1099                   DO islope_loop=islope,iflat,-1
    1100                      if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
    1101                             tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
    1102                             bool_sublim=1
    1103                             exit
    1104                      endif
    1105                   enddo
    1106                   if (bool_sublim.eq.1) then
    1107                    exit
     931          if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-5) THEN !co2ice disappeared, look for closest point without co2ice
     932            if(latitude_deg(ig).gt.0) then
     933              do ig_loop=ig,ngrid
     934                DO islope_loop=islope,iflat,-1
     935                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
     936                    tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
     937                    bool_sublim=1
     938                    exit
    1108939                  endif
    1109940                enddo
    1110               else
    1111                do ig_loop=ig,1,-1
    1112                   DO islope_loop=islope,iflat
    1113                      if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
    1114                             tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
    1115                             bool_sublim=1
    1116                             exit
    1117                      endif
    1118                   enddo
    1119                   if (bool_sublim.eq.1) then
    1120                    exit
     941                if (bool_sublim.eq.1) then
     942                  exit
     943                endif
     944              enddo
     945            else
     946              do ig_loop=ig,1,-1
     947                DO islope_loop=islope,iflat
     948                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
     949                    tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
     950                    bool_sublim=1
     951                    exit
    1121952                  endif
    1122953                enddo
    1123 
     954                if (bool_sublim.eq.1) then
     955                  exit
     956                endif
     957              enddo
     958            endif
     959            initial_co2_ice_slope(ig,islope)=0
     960            if ((co2ice_slope(ig,islope).lt.1e-5).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then   
     961              albedo_slope(ig,1,islope) = albedo_h2o_frost
     962              albedo_slope(ig,2,islope) = albedo_h2o_frost
     963            else
     964              albedo_slope(ig,1,islope) = albedodat(ig)
     965              albedo_slope(ig,2,islope) = albedodat(ig)
     966              emiss_slope(ig,islope) = emissiv         
     967            endif
     968          elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2
     969            ave=0.
     970            do t=1,timelen
     971              if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then
     972                ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
     973              else
     974                ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
    1124975              endif
    1125               initial_co2_ice_slope(ig,islope)=0
    1126 
    1127               if ((co2ice_slope(ig,islope).lt.1e-5).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then
    1128                  
    1129                albedo_slope(ig,1,islope) = albedo_h2o_frost
    1130                albedo_slope(ig,2,islope) = albedo_h2o_frost
    1131 
    1132               else
    1133                            
    1134                albedo_slope(ig,1,islope) = albedodat(ig)
    1135                albedo_slope(ig,2,islope) = albedodat(ig)
    1136                emiss_slope(ig,islope) = emissiv         
    1137 
    1138               endif
    1139 
    1140            elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2
    1141  
    1142              ave=0.
    1143              do t=1,timelen
    1144               if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then
    1145                  ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
    1146               else
    1147                  ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
    1148               endif
    1149              enddo
    1150 
    1151              tsurf_ave_phys(ig,islope)=ave/timelen
    1152            endif
     976            enddo
     977            tsurf_ave_phys(ig,islope)=ave/timelen
     978          endif
    1153979        enddo
    1154        enddo
    1155 
    1156       print *, "7 is okay"
    1157    
    1158      do t = 1,timelen
    1159      tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:))
    1160      enddo
    1161 ! for the start
    1162      do ig = 1,ngrid
    1163       do islope = 1,nslope
    1164      tsurf_slope(ig,islope) =  tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
    1165      enddo
    1166      enddo
    1167 
    1168 ! II_c.2 Update soil temperature
    1169 
    1170 allocate(TI_locslope(ngrid,nsoilmx_PEM))
    1171 allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
    1172 allocate(Tsurf_locslope(ngrid))
    1173 allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
    1174 allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
    1175 
    1176 
    1177 
    1178 
    1179 print *,"8 is ok"
     980      enddo
     981
     982      do t = 1,timelen
     983        tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:))
     984      enddo
     985      ! for the start
     986      do ig = 1,ngrid
     987        do islope = 1,nslope
     988          tsurf_slope(ig,islope) =  tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
     989        enddo
     990      enddo
     991
     992    if(soil_pem) then
     993
     994! II_d.2 Update soil temperature
     995
     996  allocate(TI_locslope(ngrid,nsoilmx_PEM))
     997  allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
     998  allocate(Tsurf_locslope(ngrid))
     999  allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
     1000  allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
     1001
     1002  print *,"Updating soil temperature"
    11801003
    11811004! Soil averaged
    … …  
    11861009     alph_locslope(:,:) = alph_PEM(:,:,islope)
    11871010     beta_locslope(:,:) = beta_PEM(:,:,islope)
    1188      call soil_pem(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
    1189      call soil_pem(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
     1011     call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
     1012     call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
    11901013     do t = 1,timelen
    11911014       tsoil_phys_PEM_timeseries(:,:,islope,t) =  tsoil_phys_PEM_timeseries(:,:,islope,t) + (Tsoil_locslope(:,:)- tsoil_PEM(:,:,islope))
    … …  
    11981021  enddo
    11991022
    1200 
    1201       print *, "9 is okay"
    1202 
    1203 
    1204 !!!!!! Check for the time series
    1205 
    1206    
    1207       do ig = 1,ngrid
     1023  print *, "Update of soil temperature done"
     1024
     1025! Check Nan in the time series
     1026     do ig = 1,ngrid
    12081027       do islope = 1,nslope
    12091028         do iloop = 1,nsoilmx_PEM
    1210  
    1211          if(isnan(tsoil_PEM(ig,iloop,islope))) then
    1212           write(*,*) "is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope)
     1029           if(isnan(tsoil_PEM(ig,iloop,islope))) then
     1030            write(*,*) "Problem : There is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope)
    12131031           stop
    1214          endif
     1032           endif
    12151033         enddo
    1216         enddo
    1217       enddo
    1218 
    1219 
    1220 write(*,*) "before ice table, no stop"
    1221 
    1222 ! II_c.3 Update the ice table
     1034       enddo
     1035     enddo
     1036
     1037  write(*,*) "Compute ice table"
     1038
     1039! II_d.3 Update the ice table
    12231040     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, &
    12241041                           ps_phys_timeseries, ice_depth)
    12251042       
    1226       print *, "10 is okay"
    1227 ! II_c.3 Update the soil thermal properties
    1228       call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),ps_phys, &
    1229         cell_area,ice_depth,TI_PEM)
    1230 
    1231 ! II_c.4 Update the mass of the regolith adsorbded
    1232    
    1233  
     1043      print *, "Update soil propreties"
     1044! II_d.4 Update the soil thermal properties
     1045      call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
     1046        ice_depth,TI_PEM)
     1047
     1048! II_d.5 Update the mass of the regolith adsorbded
    12341049   call regolith_co2adsorption(ngrid,nslope,nsoilmx_PEM,timelen,ps_phys_timeseries,tsoil_PEM,TI_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope, &
    12351050                               co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), q_co2_PEM_phys,q_h2o_PEM_phys,co2_adsorbded_phys,delta_co2_adsorbded)
    12361051
     1052  endif !soil_pem
    12371053
    12381054!------------------------
    … …  
    12401056! II  Run
    12411057!    II_a update pressure, ice and tracers
    1242 !    II_b CO2 glaciers flows
    1243 !    II_c Update surface and soil temperatures
    1244 !    II_d Update the tendencies and stopping criterion
    1245 
    1246 !------------------------
    1247 
    1248 
     1058!    II_b Evolution of the ice
     1059!    II_c CO2 glaciers flows
     1060!    II_d Update surface and soil temperatures
     1061!    II_e Update the tendencies
     1062
     1063!------------------------
     1064
     1065      print *, "Adaptation of the new co2 tendencies to the current pressure"
    12491066      call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,&     
    12501067                               global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
    12511068
    1252           print *, "12 is okay"
     1069!------------------------
     1070
     1071! II  Run
     1072!    II_a update pressure, ice and tracers
     1073!    II_b Evolution of the ice
     1074!    II_c CO2 glaciers flows
     1075!    II_d Update surface and soil temperatures
     1076!    II_e Update the tendencies
     1077!    II_f Checking the stopping criterion
     1078
     1079!------------------------
     1080      call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
     1081
     1082      call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
     1083
    12531084      year_iter=year_iter+dt_pem
    12541085
    1255 !We check with we should stop
    1256 
    1257 
    1258       call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
    1259 
    1260 
    1261           print *, "13 is okay"
    1262       call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
    1263 
    1264  print *, "criterion ok"
    1265 
    1266 
    1267       print *, "STOPPING_water", STOPPING_water, "STOPPING1_water", STOPPING_1_water
    1268       print *, "STOPPING_co2", STOPPING_co2, "STOPPING1_co2", STOPPING_1_co2
    1269       print *, "year_iter=", year_iter
     1086      print *, "Checking all the stopping criterion."
     1087      if (STOPPING_water) then
     1088        print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
     1089      endif
     1090      if (STOPPING_1_water) then
     1091        print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
     1092      endif
     1093      if (STOPPING_co2) then
     1094        print *, "STOPPING because surface of co2 ice sublimating is too low or global pressure changed too much, see message above", STOPPING_co2
     1095      endif
     1096      if (STOPPING_1_co2) then
     1097        print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
     1098      endif
     1099
    12701100      if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2)  then
    12711101        exit
     1102      else
     1103        print *, "We continue!"
     1104        print *, "Number of iteration of the PEM : year_iter=", year_iter 
    12721105      endif
    12731106
    12741107      global_ave_press_old=global_ave_press_new
    1275       zplev_old=zplev_new
    1276 
    1277 
    1278       enddo  ! big loop of the pem
     1108
     1109      enddo  ! big time iteration loop of the pem
    12791110
    12801111
    12811112!---------------------------- END RUN PEM ---------------------
    12821113
    1283 
    1284 
    1285 
    12861114!---------------------------- OUTPUT ---------------------
    12871115
    1288 
    1289 
    1290 
    12911116!------------------------
    12921117
    12931118! III Output
    1294 !    III_a Recomp tendencies for the start and startfi
    1295 
    1296 !------------------------
    1297 
    1298 
    1299 ! III_a.1 Ice update
     1119!    III_a Update surface value for the PCM start files
     1120
     1121!------------------------
     1122
     1123! III_a.1 Ice update (for startfi)
     1124! Co2 ice
    13001125      DO ig = 1,ngrid
    1301           co2ice(ig) = 0.
    1302              DO islope = 1,nslope
    1303                  co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) &
    1304                             * subslope_dist(ig,islope) /          &
    1305                            cos(pi*def_slope_mean(islope)/180.)
    1306              ENDDO
     1126        co2ice(ig) = 0.
     1127        DO islope = 1,nslope
     1128          co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) &
     1129                      * subslope_dist(ig,islope) /          &
     1130                      cos(pi*def_slope_mean(islope)/180.)
     1131        ENDDO
    13071132      ENDDO ! of DO ig=1,ngrid
    1308 
    1309       DO ig = 1,ngrid !!!!! TO DOC GHANGE IF QSURF ?ü!?!?!?
    1310           qsurf(ig,igcm_h2o_ice) = 0.
    1311              DO islope = 1,nslope
    1312                  qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
    1313                             * subslope_dist(ig,islope) /          &
    1314                            cos(pi*def_slope_mean(islope)/180.)
    1315              ENDDO
     1133! H2o ice
     1134      DO ig = 1,ngrid
     1135        qsurf(ig,igcm_h2o_ice) = 0.
     1136        DO islope = 1,nslope
     1137          qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
     1138                                   * subslope_dist(ig,islope) /          &
     1139                                  cos(pi*def_slope_mean(islope)/180.)
     1140        ENDDO
    13161141      ENDDO ! of DO ig=1,ngrid
    13171142
    1318 
    1319 
    1320 
    1321 ! III_a.2 Tsurf and Tsoil update
    1322 
    1323 
    1324    call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
    1325    tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
    1326 
     1143! III_a.2 Tsurf and Tsoil update (for startfi)
     1144
     1145   if(soil_pem) then
     1146     call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
     1147     tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
     1148   else
     1149      TI_GCM_phys(:,:,:)=TI_GCM_start(:,:,:)
     1150   endif !soil_pem
    13271151
    13281152      DO ig = 1,ngrid
    1329           tsurf(ig) = 0.
    1330          
    1331              DO islope = 1,nslope
    1332                  tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) &
    1333                             * subslope_dist(ig,islope)     
    1334              ENDDO
     1153        tsurf(ig) = 0.
     1154        DO islope = 1,nslope
     1155          tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) &
     1156                     * subslope_dist(ig,islope)     
     1157        ENDDO
    13351158      ENDDO ! of DO ig=1,ngrid
    13361159
    13371160      DO ig = 1,ngrid
    1338          DO iloop = 1,nsoilmx
     1161        DO iloop = 1,nsoilmx
    13391162          tsoil(ig,iloop) = 0.
    13401163          inertiesoil(ig,iloop) = 0.
    1341              DO islope = 1,nslope
    1342                   tsoil(ig,iloop) =  tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
     1164          DO islope = 1,nslope
     1165            tsoil(ig,iloop) =  tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
    13431166                            * subslope_dist(ig,islope)   
    1344                   inertiesoil(ig,iloop) =  inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
     1167            inertiesoil(ig,iloop) =  inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
    13451168                            * subslope_dist(ig,islope)         
    1346              ENDDO
    13471169          ENDDO
     1170        ENDDO
    13481171      ENDDO ! of DO ig=1,ngrid
    13491172
    1350 ! III_a.3 Surface optical properties
     1173! III_a.3 Surface optical properties (for startfi)
    13511174
    13521175    DO ig = 1,ngrid
    … …  
    13571180                         *subslope_dist(ig,islope)
    13581181        ENDDO
    1359        ENDDO
    13601182      ENDDO
    1361 
     1183    ENDDO
    13621184
    13631185    DO ig = 1,ngrid
    1364 
    1365         emis(ig) =0.
    1366         DO islope = 1,nslope
    1367           emis(ig)= emis(ig)+emiss_slope(ig,islope) &
     1186      emis(ig) =0.
     1187      DO islope = 1,nslope
     1188        emis(ig)= emis(ig)+emiss_slope(ig,islope) &
    13681189                         *subslope_dist(ig,islope)
    1369         ENDDO
    1370        ENDDO
    1371      
    1372 
    1373 
    1374 !------------------------
     1190      ENDDO
     1191    ENDDO
     1192
     1193! III_a.4 Pressure (for start)
     1194     do i=1,ip1jmp1
     1195       ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM
     1196     enddo
     1197
     1198     do i = 1,ngrid
     1199       ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_GCM
     1200     enddo
     1201
     1202! III_a.5 Tracer (for start)
     1203     allocate(zplev_new(ngrid,nlayer+1))
     1204
     1205     do l=1,nlayer+1
     1206       do ig=1,ngrid
     1207         zplev_new(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
     1208       enddo
     1209     enddo
     1210
     1211     DO nnq=1,nqtot
     1212       if (noms(nnq).NE."co2") then
     1213         DO l=1,llm-1
     1214           DO ig=1,ngrid
     1215             q(ig,l,nnq)=q(ig,l,nnq)*(zplev_gcm(ig,l)-zplev_gcm(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))
     1216           ENDDO
     1217           q(:,llm,nnq)=q(:,llm-1,nnq)
     1218         ENDDO
     1219       else
     1220         DO l=1,llm-1
     1221           DO ig=1,ngrid
     1222             q(ig,l,nnq)=q(ig,l,nnq)*(zplev_gcm(ig,l)-zplev_gcm(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))  &
     1223                             + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
     1224                                   (zplev_gcm(ig,l)-zplev_gcm(ig,l+1))     )  / &
     1225                                   (zplev_new(ig,l)-zplev_new(ig,l+1))
     1226           ENDDO
     1227           q(:,llm,nnq)=q(:,llm-1,nnq)
     1228         ENDDO
     1229       endif
     1230     ENDDO
     1231
     1232! Conserving the tracers's mass for GCM start files
     1233     DO nnq=1,nqtot
     1234       DO ig=1,ngrid
     1235         DO l=1,llm-1
     1236           if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number" .or. noms(nnq).NE."ccn_number") ) then
     1237             extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
     1238             q(ig,l,nnq)=1.
     1239             q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
     1240           endif
     1241           if(q(ig,l,nnq).LT.0) then
     1242             q(ig,l,nnq)=1E-30
     1243           endif
     1244         ENDDO
     1245       enddo
     1246     enddo
     1247     
     1248!------------------------
     1249     if(evol_orbit_pem) then
     1250      call recomp_orb_param(year_iter)
     1251     endif
    13751252
    13761253! III Output
    1377 !    III_a Recomp tendencies for the start and startfi
     1254!    III_a Update surface value for the PCM start files
    13781255!    III_b Write start and starfi.nc
    13791256
    … …  
    13861263      ztime_fin=0.
    13871264
    1388          print *, "RVip1jmp1", ip1jmp1
    13891265     allocate(p(ip1jmp1,nlayer+1))
    1390          print *, "ngrid", ngrid
    1391          print *, "nlayer", nlayer
    1392           CALL pression (ip1jmp1,ap,bp,ps,p)
    1393      print *, "masse 1", masse(1,:)
    1394           CALL massdair(p,masse)
    1395      print *, "masse 2", masse(1,:)
    1396 
    1397      do nnq=1,nqtot
    1398         call gr_fi_dyn(llm,ngridmx,iip1,jjp1,q_phys(:,:,nnq),q(:,:,nnq))
    1399      enddo
    1400 
     1266     CALL pression (ip1jmp1,ap,bp,ps,p)
     1267     CALL massdair(p,masse)
    14011268
    14021269      CALL dynredem0("restart_evol.nc", day_ini, phis)
    … …  
    14041271      CALL dynredem1("restart_evol.nc",   &
    14051272                time_0,vcov,ucov,teta,q,masse,ps)
    1406 
     1273      print *, "restart_evol.nc has been written"
    14071274
    14081275! III_b.2 WRITE restartfi.nc
    1409 
    14101276     
    1411              call physdem0("restartfi_evol.nc",longitude,latitude, &
     1277      call physdem0("restartfi_evol.nc",longitude,latitude, &
    14121278                        nsoilmx,ngrid,nlayer,nq,              &
    14131279                        ptimestep,pday,0.,cell_area,          &
    … …  
    14161282                        subslope_dist)
    14171283
    1418 
    1419           call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
     1284     call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
    14201285                     ptimestep,ztime_fin,                        &
    14211286                     tsurf,tsoil,co2ice,albedo,emis,             &
    … …  
    14251290                     emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys)
    14261291
    1427 
     1292      print *, "restartfi_evol.nc has been written"
    14281293!------------------------
    14291294
    14301295! III Output
    1431 !    III_a Recomp tendencies for the start and startfi
     1296!    III_a Update surface value for the PCM start files
    14321297!    III_b Write start and starfi.nc
    14331298!    III_c Write start_pem
    … …  
    14351300!------------------------
    14361301
    1437  
    1438 
    1439          call pemdem1("restartfi_PEM.nc",ztime_fin,nsoilmx_PEM,ngrid,nslope , &
     1302         call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
    14401303                      tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys)
    14411304
    1442 
    1443  
    1444 
     1305      print *, "restartfi_PEM.nc has been written"
     1306      print *, "The PEM had run for ", year_iter, " years."
    14451307      print *, "LL & RV : So far so good"
    14461308
    1447 
    1448 
    1449 
    1450 
    1451 
    14521309END PROGRAM pem
    1453 

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

    r2794 r2835  
    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_yr1_inst,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, m_co2_regolith_phys,deltam_co2_regolith_phys)
    33   
    4    use iostart_PEM, only:  open_startphy, close_startphy, get_field 
     4   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
    55   use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
    66   use comsoil_h, only:  volcapa,inertiedat
    77   use ini_soil_mod, only:  ini_icetable
    8    use soil_evolution_mod, only: soil_pem_CN
     8   use soil_evolution_mod, only: soil_pem,soil_pem_CN
    99   use adsorption_mod, only : regolith_co2adsorption
     10   USE temps_mod_evol, ONLY: year_PEM
     11
    1012  implicit none
    11 
    1213 
    1314  character(len=*), intent(in) :: filename
    … …  
    2324  real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
    2425  real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
    25   real,intent(in) :: ps_yr1_inst(ngrid,timelen)
    2626  real,intent(in) :: ps_inst(ngrid,timelen)                        ! surface pressure [Pa]
    2727  real,intent(in) :: tsurf_inst(ngrid,nslope,timelen) ! soil (mid-layer) temperature
    … …  
    3434
    3535! outputs
    36 
    37 
    38 
    39 
    4036
    4137  real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) thermal inertia (SI)
    … …  
    5955   real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)
    6056   real :: tsoil_tmp(ngrid,nsoil_PEM,nslope)
    61    real :: alph_tmp(nsoil_PEM-1)
    62    real :: beta_tmp(nsoil_PEM-1)
    63    real :: co2_ads_prev(ngrid) 
    64    real :: ps_ave_yr1(ngrid)
    65    real :: ps_ave_yr2(ngrid)
     57   real :: alph_tmp(ngrid,nsoil_PEM-1)
     58   real :: beta_tmp(ngrid,nsoil_PEM-1)
     59   real :: co2_ads_prev(ngrid)
     60   real :: year_PEM_read
     61
    6662!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    6763!!!
    … …  
    7874!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    7975
    80 
    81 ! 0. Some initializations
    82  ps_ave_yr1(:) = sum(ps_yr1_inst(:,:),2)/timelen
    83  ps_ave_yr2(:) = sum(ps_inst(:,:),2)/timelen
    84 
    8576!1. Run
    86 
    87 
    88 
    8977
    9078if (startpem_file) then
    9179   ! open pem initial state file:
    9280   call open_startphy(filename)
     81
     82    if(soil_pem) then
     83
     84   call get_var("Time",year_PEM_read,found)
     85   year_PEM=INT(year_PEM_read)
     86   if(.not.found) then
     87      write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
     88   else
     89      write(*,*)'year_PEM of startpem=', year_PEM
     90   endif
    9391 
    9492!1. Thermal Inertia
    … …  
    103101
    104102      do ig = 1,ngrid
    105 
    106103          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then
    107104!!! transition
    … …  
    110107            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
    111108                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
    112 
    113           do iloop=nsoil_GCM+2,n_1km
    114             TI_PEM(ig,iloop,islope) = TI_breccia
    115          enddo     
    116          else ! we keep the high ti values
    117            do iloop=nsoil_GCM+1,n_1km
    118                   TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
    119            enddo
    120 
    121          endif ! TI PEM and breccia comparison
    122 
    123 
     109            do iloop=nsoil_GCM+2,n_1km
     110              TI_PEM(ig,iloop,islope) = TI_breccia
     111            enddo     
     112          else ! we keep the high ti values
     113            do iloop=nsoil_GCM+1,n_1km
     114              TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
     115            enddo
     116          endif ! TI PEM and breccia comparison
    124117!! transition
    125              delta = 1000.
    126              TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
    127             (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
    128                       ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
     118          delta = 1000.
     119          TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
     120               (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
     121               ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
    129122          do iloop=n_1km+2,nsoil_PEM
    130123            TI_PEM(ig,iloop,islope) = TI_bedrock
    131          enddo   
    132       enddo
    133 
     124          enddo   
     125      enddo
    134126   else
    135127     do iloop = nsoil_GCM+1,nsoil_PEM
    … …  
    141133print *,'PEMETAT0: THERMAL INERTIA DONE'
    142134
    143 
    144 
    145135! b. Special case for inertiedat, inertiedat_PEM
    146136call get_field("inertiedat_PEM",inertiedat_PEM,found)
    … …  
    148138 do iloop = 1,nsoil_GCM
    149139   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
    150  enddo
    151 
    152  
     140 enddo
    153141!!! zone de transition
    154142delta = 50.
    … …  
    159147                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
    160148
    161 
    162149 do iloop = nsoil_GCM+2,n_1km
    163150       inertiedat_PEM(ig,iloop) = TI_breccia
    … …  
    179166enddo
    180167
    181 
    182 
    183168 do iloop = n_1km+2, nsoil_PEM
    184169   do ig = 1,ngrid
    … …  
    187172 enddo
    188173endif ! not found
    189 
    190174
    191175!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    221205
    222206    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
    223     call soil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
    224     call soil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
     207    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
     208    call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
    225209    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
    226     call soil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
     210    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
    227211
    228212
    … …  
    237221        enddo
    238222     enddo
    239 
    240 
    241 
    242223   
    243224ENDDO
     225
    244226print *,'PEMETAT0: SOIL TEMP  DONE'
    245227
    … …  
    248230!3. Ice Table
    249231
    250  
    251232   call get_field("ice_table",ice_table,found)
    252233   if(.not.found) then
    … …  
    254235      write(*,*)'will reconstruct the values of ice table'
    255236
    256 
    257237      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))
    258 
    259238
    260239   else
    … …  
    262241     call computeice_table(timelen,ngrid,nslope,nsoil_PEM,tsoil_inst,tsurf_inst,q_co2,q_h2o, ps_inst, ice_table)
    263242
    264 
    265243   endif
    266244
    267245print *,'PEMETAT0: ICE TABLE  DONE'
    268246
    269 
    270 
    271 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    272 
     247!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    273248
    274249!4. CO2 Adsorption
    … …  
    299274print *,'PEMETAT0: CO2 adsorption done  '
    300275
    301 
     276endif ! soil_pem
    302277
    303278call close_startphy
    304279
    305 
    306280else !No startfi, let's build all by hand
    307281
    308 
     282    year_PEM=0
     283
     284    if(soil_pem) then
    309285
    310286!a) Thermal inertia
    … …  
    326302                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
    327303           enddo
    328 
    329304         endif
    330 
    331305
    332306!! transition
    … …  
    339313         enddo   
    340314      enddo
    341 
    342315enddo
    343 
    344316
    345317 do iloop = 1,nsoil_GCM
    … …  
    440412print *,'PEMETAT0: CO2 adsorption done  '
    441413
    442 
    443 
    444 
    445 
    446 
    447 
    448 
    449 
     414endif !soil_pem
    450415
    451416endif ! of if (startphy_file)
    452417
    453418!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     419   
     420     if(soil_pem) then
     421
    454422    DO ig = 1,ngrid
    455423      DO islope = 1,nslope
    … …  
    459427      ENDDO
    460428    ENDDO
    461 
    462 
    463429
    464430!! small test
    … …  
    474440     ENDDO
    475441
     442     endif!soil_pem
     443
    476444      write(*,*) "construction ok, no nan"
    477445     

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

    r2794 r2835  
    88                         day_ini,time,nslope,def_slope,subslope_dist)
    99! create physics restart file and write time-independent variables
    10 !  use comsoil_h, only: inertiedat, volcapa, mlayer
    1110  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM
    1211  use iostart_PEM, only : open_restartphy, close_restartphy, &
    … …  
    3029  real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes
    3130  real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics
    32 
    33 
    3431 
    3532  ! Create physics start file
    … …  
    5956end subroutine pemdem0
    6057
    61 subroutine pemdem1(filename,time,nsoil_PEM,ngrid,nslope, &
     58subroutine pemdem1(filename,year_iter,nsoil_PEM,ngrid,nslope, &
    6259                    tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table,m_co2_regolith)
    6360  ! write time-dependent variable to restart file
    … …  
    6663
    6764  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM
     65  USE temps_mod_evol, ONLY: year_PEM
     66  use soil_evolution_mod, only: soil_pem
    6867               
    6968  implicit none
    … …  
    7473  integer,intent(in) :: nsoil_PEM
    7574  integer,intent(in) :: ngrid
    76   real,intent(in) :: time
    7775  real,intent(IN) :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
    7876  real,intent(IN) :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
    … …  
    8280  integer :: islope
    8381  CHARACTER*2 :: num 
     82  integer, intent(IN) :: year_iter
    8483
     84  real :: year_tot
    8585 
    8686  integer :: iq
    … …  
    8888  ! indexes of water vapour & water ice tracers (if any):
    8989
    90  
     90  year_tot=real(year_iter+year_PEM)
     91  print *, "Year of simulation=", year_tot 
    9192
    9293  ! Open file
    … …  
    9596  ! First variable to write must be "Time", in order to correctly
    9697  ! set time counter in file
    97   call put_var("Time","Temps de simulation",time)
    98  
    99  
     98  call put_var("Time","Year of simulation",year_tot)
     99
     100    if(soil_pem) then
    100101 
    101102  ! Multidimensionnal variables (undermesh slope statistics)
    … …  
    103104  write(num,fmt='(i2.2)') islope
    104105  call put_field("tsoil_PEM_slope"//num,"Soil temperature by slope type", &
    105                  tsoil_slope_PEM(:,:,islope),time)
     106                 tsoil_slope_PEM(:,:,islope),year_tot)
    106107  call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type", &
    107                  inertiesoil_slope_PEM(:,:,islope),time)
     108                 inertiesoil_slope_PEM(:,:,islope),year_tot)
    108109
    109110  call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
    110                     m_co2_regolith(:,:,islope), time)
     111                    m_co2_regolith(:,:,islope), year_tot)
    111112
    112113ENDDO
    113114
    114115  call put_field("ice_table","Depth of ice table", &
    115                  ice_table,time)
     116                 ice_table,year_tot)
    116117
    117118  call put_field("inertiedat_PEM","Thermal inertie of PEM ", &
    118                  inertiedat_PEM,time)
     119                 inertiedat_PEM,year_tot)
    119120
     121   endif !soil_pem
    120122
    121123  ! Close file

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

    r2794 r2835  
    99                        nf90_inquire_dimension,nf90_close
    1010      use comsoil_h, only: nsoilmx
     11      USE soil_evolution_mod, ONLY: soil_pem
    1112
    1213      IMPLICIT NONE
    … …  
    2324! Arguments:
    2425  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
     26  INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
    2527
    2628  INTEGER :: iim_input,jjm_input,nlayer,nslope
    … …  
    3436  REAL, INTENT(OUT) ::  min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
    3537  REAL, INTENT(OUT) ::  min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
    36 !  REAL, ALLOCATABLE ::  vmr_co2_gcm(:,:,:)                     !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
    37   REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,2676)                     !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
    38 !  REAL, ALLOCATABLE ::  q_h2o_GCM(:,:,:)
    39   REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,2676)
    40   REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,2676)
    41 !  REAL, ALLOCATABLE ::  q_co2_GCM(:,:,:)
     38  REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)      !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
     39  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,timelen)
     40  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,timelen)
    4241  REAL, ALLOCATABLE ::  q1_co2_GCM(:,:,:)
    43 !  real, INTENT(OUT) ::  vmr_co2_phys_gcm(:,:)                  !!!!vmr_co2_gcm(ngrid,timelen)
    44 !  REAL, ALLOCATABLE ::  ps_GCM(:,:,:)
    45   REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,2676)
    46 
     42  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,timelen)
    4743
    4844!SOIL
    … …  
    5046  REAL, INTENT(OUT) ::  tsoil_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average soil temperature of the concatenated file
    5147
    52   REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,2676) ! Surface temperature of the concatenated file
    53   REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,2676) ! Soil temperature of the concatenated file
    54 
    55   REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,2676) ! Thermal Inertia  of the concatenated file
     48  REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file
     49  REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
     50
     51  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Thermal Inertia  of the concatenated file
    5652  REAL, INTENT(OUT) ::  TI_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average Thermal Inertia  of the concatenated file
    57   REAL, INTENT(OUT) ::  co2_ice_slope(iim_input+1,jjm_input+1,nslope,2676) ! Minimum of co2_ice slope of the year
     53  REAL, INTENT(OUT) ::  co2_ice_slope(iim_input+1,jjm_input+1,nslope,timelen) ! Minimum of co2_ice slope of the year
    5854!===============================================================================
    5955!   Local Variables
    … …  
    6561
    6662  REAL,ALLOCATABLE :: time(:) ! times stored in start
    67   INTEGER :: timelen ! number of times stored in the file
    6863  INTEGER :: indextime ! index of selected time
    6964
    … …  
    7469  INTEGER :: islope
    7570  CHARACTER*2 :: num
    76 
    7771
    7872!-----------------------------------------------------------------------
    … …  
    8478      B=1/m_noco2
    8579
     80      allocate(co2_ice_s(iim+1,jjm+1,timelen))
     81      allocate(q1_co2_GCM(iim+1,jjm+1,timelen))
     82      allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
     83      allocate(h2o_ice_s(iim+1,jjm+1,timelen))
     84
     85  print *, "Opening ", fichnom, "..."
     86
    8687!  Open initial state NetCDF file
    8788  var=fichnom
    8889  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
    8990
    90       ierr = nf90_inq_varid (fID, "temps", vID)
    91       IF (ierr .NE. nf90_noerr) THEN
    92         write(*,*)"pemetat0: Le champ <temps> est absent"
    93         write(*,*)"pemetat0: J essaie <Time>"
    94         ierr = nf90_inq_varid (fID, "Time", vID)
    95         IF (ierr .NE. nf90_noerr) THEN
    96            write(*,*)"pemetat0: Le champ <Time> est absent"
    97            write(*,*)trim(nf90_strerror(ierr))
    98      print *, "ICIIII0"
    99            CALL ABORT_gcm("pemetat0", "", 1)
    100         ENDIF
    101         ! Get the length of the "Time" dimension
    102      print *, "ICIIIITIME"
    103         ierr = nf90_inq_dimid(fID,"Time",vID)
    104         ierr = nf90_inquire_dimension(fID,vID,len=timelen)
    105       ELSE   
    106         ! Get the length of the "temps" dimension
    107      print *, "ICIIIITEMPS"
    108         ierr = nf90_inq_dimid(fID,"temps",vID)
    109         ierr = nf90_inquire_dimension(fID,vID,len=timelen)
    110       ENDIF
    111 
    112      print *, "ICIIII"
    113 
    114       allocate(co2_ice_s(iim+1,jjm+1,timelen))
    115 
    116      print *, "ICIIIIAAAA"
    117 
    118 !      allocate(q_co2_GCM(iim+1,jjm+1,timelen))
    119 
    120      print *, "ICIIIIBBBB"
    121 
    122 !      allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
    123 
    124      print *, "ICIIIICCCC"
    125 
    126       allocate(q1_co2_GCM(iim+1,jjm+1,timelen))
    127 
    128      print *, "ICIIII2"
    129 
    130 
    131 
    132       allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
    133       allocate(h2o_ice_s(iim+1,jjm+1,timelen))
    134           print *, "ICIIII3"
     91     print *, "Downloading data for h2oice ..."
    13592
    13693! Get h2o_ice_s of the concatenated file
    13794  CALL get_var3("h2o_ice_s"   ,h2o_ice_s)
    13895
    139   print *, "A"
     96     print *, "Downloading data for h2oice done"
     97     print *, "Downloading data for co2ice ..."
    14098
    14199  CALL get_var3("co2ice"   ,co2_ice_s)
     100
     101     print *, "Downloading data for co2ice done"
     102     print *, "Downloading data for vmr co2..."
     103
    142104  CALL get_var3("co2_cropped"   ,q_co2_GCM)
     105
     106     print *, "Downloading data for vmr co2 done"
     107     print *, "Downloading data for vmr h20..."
     108
    143109  CALL get_var3("h2o_cropped"   ,q_h2o_GCM)
    144110
    145   print *, "B"
     111     print *, "Downloading data for vmr h2o done"
     112     print *, "Downloading data for surface pressure ..."
    146113
    147114  CALL get_var3("ps"   ,ps_GCM)
    148115
    149   print *, "C"
    150 
    151   print *, "nslope", nslope
     116     print *, "Downloading data for surface pressure done"
     117     print *, "nslope=", nslope
     118     print *, "Downloading data for co2ice_slope ..."
    152119
    153120DO islope=1,nslope
    … …  
    156123ENDDO
    157124
    158   print *, "co2ice_slope"
     125     print *, "Downloading data for co2ice_slope done"
     126     print *, "Downloading data for h2o_ice_s_slope ..."
    159127
    160128DO islope=1,nslope
    … …  
    163131ENDDO
    164132
    165   print *, "h2o_ice_s_slope"
     133     print *, "Downloading data for h2o_ice_s_slope done"
     134     print *, "Downloading data for tsurf_slope ..."
    166135
    167136DO islope=1,nslope
    … …  
    170139ENDDO
    171140
    172   print *, "tsurf_slope"
     141     print *, "Downloading data for tsurf_slope done"
     142
     143     if(soil_pem) then
     144
     145     print *, "Downloading data for tsoil_slope ..."
    173146
    174147DO islope=1,nslope
    … …  
    177150ENDDO
    178151
    179   print *, "tsoil_slope"
     152     print *, "Downloading data for tsoil_slope done"
     153     print *, "Downloading data for inertiesoil_slope ..."
    180154
    181155DO islope=1,nslope
    … …  
    184158ENDDO
    185159
    186   print *, "inertiesoil_slope"
    187 
    188 
    189 
    190 
     160     print *, "Downloading data for inertiesoil_slope done"
     161
     162  endif
    191163
    192164! Compute the minimum over the year for each point
     165  print *, "Computing the min of h2o_ice"
    193166  min_h2o_ice_s(:,:)=minval(h2o_ice_s,3)
     167  print *, "Computing the min of co2_ice"
    194168  min_co2_ice_s(:,:)=minval(co2_ice_s,3)
    195169
     170  print *, "Computing the min of h2o_ice_slope"
     171  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
     172  print *, "Computing the min of co2_ice_slope"
    196173  min_co2_ice_slope(:,:,:)=minval(co2_ice_slope,4)
    197   min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
    198174
    199175!Compute averages
    200176
    201 !  DO i=1,timelen
     177    print *, "Computing average of tsurf"
    202178    tsurf_ave(:,:,:)=SUM(tsurf_gcm(:,:,:,:),4)/timelen
     179
     180  if(soil_pem) then
     181    print *, "Computing average of tsoil"
    203182    tsoil_ave(:,:,:,:)=SUM(tsoil_gcm(:,:,:,:,:),5)/timelen
     183    print *, "Computing average of TI"
    204184    TI_ave(:,:,:,:)=SUM(TI_gcm(:,:,:,:,:),5)/timelen
    205 !  ENDDO
    206 
    207 
    208 
     185  endif
    209186
    210187! By definition, a density is positive, we get rid of the negative values
    … …  
    229206    DO j = 1, jjm+1
    230207      DO t = 1, timelen
    231 !         q1_co2_GCM(i,j,t)=q_co2_GCM(i,j,t)
    232 !         ps_GCM(i,j,t)=ps(i)
    233 !c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
    234208         if (q_co2_GCM(i,j,t).LT.0) then
    235209              q_co2_GCM(i,j,t)=1E-10
    … …  
    248222  ENDDO
    249223
    250 
    251 
    252 
    253 
    254 
    255224  CONTAINS
    256225

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

    r2794 r2835  
    2020  REAL, INTENT(in) ::  vmr_co2_pem(ngrid,timelen)                ! physical point field : Volume mixing ratio of co2 in the first layer
    2121  REAL, intent(in) :: ps_GCM_2(ngrid,timelen)                 ! physical point field : Surface pressure in the GCM
    22 !  REAL, INTENT(in) ::  q_co2_GCM(ngrid)                ! physical point field : Density of co2 in the first layer
    23 !  REAL, intent(in) :: ps_GCM(ngrid)                 ! physical point field : Density of co2 in the first layer
    2422  REAL, intent(in) :: global_ave_press_GCM
    2523  REAL, intent(in) :: global_ave_press_new
    2624  REAL, intent(in) ::  tendencies_co2_ice_phys_ini(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    2725
    28 
    2926!   OUTPUT
    3027  REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    31 
    3228
    3329!   local:
    … …  
    4541  coef=669*24*3600*eps*sigma/L
    4642
    47   print *, "coef", coef
    48   print *, "global_ave_press_GCM", global_ave_press_GCM
    49   print *, "global_ave_press_new", global_ave_press_new
    50 
    5143! Evolution of the water ice for each physical point
    5244  do i=1,ngrid
    5345    do islope=1,nslope
    54        ave=0.
    55     do t=1,timelen
    56 
    57 !       write(*,*)'i,t=',i,t,islope, alpha,beta,ave,vmr_co2_gcm(i,t),vmr_co2_pem(i,t),ps_GCM_2(i,t),global_ave_press_GCM,global_ave_press_new
    58        ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100)))**4  &
     46      ave=0.
     47      do t=1,timelen
     48        ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100)))**4  &
    5949              -(beta/(alpha-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*global_ave_press_GCM/global_ave_press_new/100)))**4
     50      enddo
     51      tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)-coef*ave/timelen
    6052    enddo
    61 !    print *, "i", i
    62 !    print *, "tendencies_co2_ice_phys_ini bef", tendencies_co2_ice_phys_ini(i,islope)
    63 !    tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)+coef*ave/timelen
    64     tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)-coef*ave/timelen
    65 
    66 !    print *, "tendencies after", tendencies_co2_ice_phys(i,islope)
    67 !    print *, "ave", ave
    68 !    print *, "timelen", timelen
    69 !    print *, "vmr_co2_pem(i,t)*ps_GCM_2(i,t)", vmr_co2_pem(i,t)*ps_GCM_2(i,t)
    70 !    print *, "-------------------"
    71   enddo
    7253  enddo
    7354

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

    r2794 r2835  
    11      SUBROUTINE soil_TIfeedback_PEM(ngrid,nsoil,icecover,   newtherm_i)
    22
    3 !      use tracer_mod, only:  rho_ice
    43      use comsoil_h_PEM, only: layer_PEM, inertiedat_PEM
    54
    … …  
    2423!  Author: Adapted from J.-B. Madeleine Mars 2008 ( Updated November 2012) by LL, 2022
    2524!=======================================================================
    26  
    2725
    2826!Local variables
    … …  
    4240
    4341      REAL ,INTENT(IN):: icecover(ngrid)         ! tracer on the surface (kg.m-2)
    44                                        
    45                                         ! water ice
     42
    4643!Outputs
    4744!-------
    4845
    4946      REAL,INTENT(INOUT) :: newtherm_i(ngrid,nsoil)    ! New soil thermal inertia
    50 
    5147
    5248      prev_thermi(:,:) = newtherm_i(:,:)

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

    r2794 r2835  
    11      MODULE soil_evolution_mod
    22
     3      IMPLICIT NONE
    34
    4       IMPLICIT NONE
     5      LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
    56     
    67      CONTAINS   
    7 
    88
    99     subroutine soil_pem_CN(ngrid,nsoil, &
    … …  
    2525
    2626#include "dimensions.h"
    27 !#include "dimphys.h"
    28 
    29 !#include"comsoil.h"
    30 
    3127
    3228!-----------------------------------------------------------------------
    … …  
    4440      real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
    4541
    46 
    47      
    4842! local variables:
    4943      integer ig,ik,isoil
    … …  
    5852      real :: Tcol(nsoil)
    5953
    60 
    6154do ig = 1,ngrid
    6255
    6356   Tcol(:) = tsoil(ig,:)
    6457!0. Build thermal properties of the ground
    65 
    6658
    6759  do isoil = 1,nsoil
    … …  
    8779beta_pem(1) = k_soil(1)*timestep/(rhoc(1)*layer_PEM(2)*layer_PEM(1))
    8880
    89 
    90 
    9181!1.c At the bottom: dT/dz = -Fgeo/k.
    9282beta_pem(nsoil) = timestep*k_soil(nsoil)/(2.*rhoc(nsoil)*(layer_PEM(isoil) - layer_PEM(isoil-1))**2)
    … …  
    10797  timestep/rhoc(nsoil)*fluxgeo/(layer_PEM(nsoil)-layer_PEM(nsoil-1))
    10898
    109 
    110 
    11199!2. Update by tridiagonal inversion
    112100  call tridag(d1,d2,d3,re,Tcol,nsoil)
    113  
    114 
    115 
    116       tsoil(ig,:) = Tcol(:)
     101  tsoil(ig,:) = Tcol(:)
    117102  enddo
    118103

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

    r2794 r2835  
    1       subroutine soil_pem(ngrid,nsoil,firstcall, &
     1      subroutine soil_pem_routine(ngrid,nsoil,firstcall, &
    22               therm_i,                          &
    33               timestep,tsurf,tsoil,alph_PEM,beta_PEM)
    … …  
    2020
    2121#include "dimensions.h"
    22 !#include "dimphys.h"
    23 
    24 !#include"comsoil.h"
    25 
    2622
    2723!-----------------------------------------------------------------------
    … …  
    4036      real,intent(inout) :: alph_PEM(ngrid,nsoil-1)
    4137      real,intent(inout) :: beta_PEM(ngrid,nsoil-1)
    42 
    43 
    44 
    45 
    4638     
    4739! local variables:
    48       integer ig,ik
    49 
    50 
    51      
     40      integer ig,ik   
    5241
    5342! 0. Initialisations and preprocessing step
    … …  
    5746      do ig=1,ngrid
    5847        do ik=0,nsoil-1
    59           mthermdiff_PEM(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
    60    
     48          mthermdiff_PEM(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa   
    6149        enddo
    6250      enddo
    63 
    6451
    6552! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
    … …  
    6956                     +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ig,ik-1))  &
    7057                         /(mlayer_PEM(ik)-mlayer_PEM(ik-1))
    71 
    7258        enddo
    7359      enddo
    … …  
    10086        enddo
    10187
    102 
    10388      enddo ! of do ig=1,ngrid
    10489
    105      
    106  
    10790endif ! of if (firstcall)
    108 
    109 
    11091
    11192      IF (.not.firstcall) THEN
    11293!  2. Compute soil temperatures
    11394! First layer:
    114 
    115 
    116 
    117 
    118       do ig=1,ngrid
    119         tsoil(ig,1)=(tsurf(ig)+mu_PEM*beta_PEM(ig,1)* & 
     95        do ig=1,ngrid
     96          tsoil(ig,1)=(tsurf(ig)+mu_PEM*beta_PEM(ig,1)* & 
    12097                                  thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
    12198                   (1.+mu_PEM*(1.0-alph_PEM(ig,1))*&
    12299                    thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
    123      
    124 
    125  
    126 
    127100
    128101! Other layers:
    129       do ik=1,nsoil-1
    130 
    131           tsoil(ig,ik+1)=alph_PEM(ig,ik)*tsoil(ig,ik)+beta_PEM(ig,ik)
    132 
    133  
    134        enddo
    135      
    136      enddo
    137           ENDIF
    138 
    139 
    140 
     102          do ik=1,nsoil-1
     103            tsoil(ig,ik+1)=alph_PEM(ig,ik)*tsoil(ig,ik)+beta_PEM(ig,ik)
     104          enddo
     105        enddo
     106      ENDIF
    141107
    142108!  2. Compute beta_PEM coefficients (preprocessing for next time step)
    … …  
    157123      enddo
    158124
    159 
    160 
    161125      end
    162126

  • trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F

    r2794 r2835  
    5656      real alpha,lay1 ! coefficients for building layers
    5757      real xmin,xmax ! to display min and max of a field
    58 
    59 
    6058     
    6159!======================================================================
    … …  
    6563
    6664! 1.4 Build mlayer(), if necessary
    67 !      if (interpol) then
    6865      ! default mlayer distribution, following a power law:
    6966      !  mlayer(k)=lay1*alpha**(k-1/2)
    … …  
    7370          mlayer_PEM(iloop)=lay1*(alpha**(iloop-0.5))
    7471        enddo
    75 !      endif
    7672
    7773! 1.5 Build layer(); following the same law as mlayer()
    … …  
    8480      enddo
    8581
    86 
    87 
    88 
    8982! 2. Thermal inertia (note: it is declared in comsoil_h)
    9083! ------------------
    91 
    9284
    9385      do ig = 1,ngrid
    … …  
    10395        enddo
    10496      enddo
     97
    10598      end subroutine soil_settings_PEM

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

    r2779 r2835  
    33IMPLICIT NONE 
    44
    5   INTEGER   nyear           !     nyear   : Maximun number of year over which the PEM can interpolate
     5  INTEGER   year_bp_ini     !     year_bp_ini : Initial year of the simulation of the PEM (in evol.def)
    66  INTEGER   dt_pem          !     dt_pem  : in years, the time step used by the PEM   
    77  REAL      alpha_criterion !     alpha_criterion : percentage of change before stopping the PEM
    8 
    9 !$OMP THREADPRIVATE(nyear)       
    10 
    11 !WARNING: when adding a threadprivate variable in this module
    12 !        do not forget to add it to the copyin clause when opening an OpenMP
    13 !        parallel section. e.g. in gcm before call leapfrog_loc and/or
    14 !        possibly in iniphysiq
     8  INTEGER   year_PEM        !     year written in startfiPEM.nc
     9  INTEGER   Max_iter_pem    !     Maximal number of iteration when converging to a steady state, read in evol.def
     10  LOGICAL   evol_orbit_pem  !     True if we want to follow the orbital parameters of ob_ex_lsp.asc, read in evol.def
    1511
    1612END MODULE temps_mod_evol

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

    r2794 r2835  
    1    SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,ps_new,&
    2                           cellarea,ice_depth,TI_PEM)
    3 
     1   SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,p_avg_new,&
     2                          ice_depth,TI_PEM)
    43
    54 USE comsoil_h, only:  inertiedat, volcapa
    … …  
    109 INTEGER,INTENT(IN) ::  ngrid, nslope, nsoil_PEM
    1110 REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope)
    12  REAL,INTENT(IN) :: ps_new(ngrid)
    13  REAL,INTENT(IN) :: cellarea(ngrid)
     11 REAL,INTENT(IN) :: p_avg_new
    1412 REAL,INTENT(IN) :: co2ice(ngrid,nslope)
    1513 REAL,INTENT(IN) :: waterice(ngrid,nslope)
    1614 REAL,INTENT(in) :: ice_depth(ngrid,nslope)
    17  
    1815
    1916! Outputs:
    … …  
    4441 REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith
    4542 REAL :: d(ngrid,nsoil_PEM,nslope)
    46  REAL :: p_avg_new
    4743 REAL :: Total_surface
    4844 INTEGER :: ispermanent_co2glaciers(ngrid,nslope)
    4945 INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)
    5046
    51 
    5247! 0. Initialisation
    53 
    54 
    55 
    56 
    57   Total_surface = sum(cellarea)
    58   p_avg_new = 0.
    59   do ig = 1,ngrid
    60     p_avg_new = p_avg_new + ps_new(ig)*cellarea(ig)/Total_surface
    61   enddo
    62 
    63 
    6448
    6549  do ig = 1,ngrid
    … …  
    8064  enddo
    8165
    82 
    8366 ispermanent_co2glaciers(:,:) = 0
    8467 ispermanent_h2oglaciers(:,:) = 0
    8568
    86 
    8769!  1.Ice TI feedback
    8870
    89 !  do ig=1,ngrid
    90 !    do islope=1,nslope
    91 !      if ((ispermanent_co2glaciers(ig,islope).eq.1).or.(ispermanent_h2oglaciers(ig,islope).eq.1)) then
    92 !       do iloop = 1,n_1km
    93 !         TI_PEM(ig,iloop,islope) = surfaceice_inertia
    94 !       enddo
    95 !      endif
    96 !    enddo
    97 !   enddo
    9871    do islope = 1,nslope
    9972     call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope),   TI_PEM(:,:,islope))
    … …  
    10174
    10275! 2. Modification of the regolith thermal inertia.
    103 
    104 
    10576
    10677  do ig=1,ngrid
    … …  
    11586   enddo
    11687enddo
    117 
    118 
    119 
    120 
    121 
    12288
    12389!  3. Build new TI for the PEM
    … …  
    142108         exit
    143109        endif
    144        
    145110      enddo
    146 
    147 
    148111
    149112      ! 4.2 Build the new ti
    … …  
    173136    enddo !ig
    174137
    175 
    176 
    177 
    178 
    179 
    180 
    181138!=======================================================================
    182139      RETURN
Note: See TracChangeset for help on using the changeset viewer.