Ignore:
Timestamp:
Sep 9, 2022, 4:02:51 PM (2 years ago)
Author:
llange
Message:

MARS PEM:

  • Add a PEMETAT0 that read "startfi_pem.nc"
  • Add the soil in the model: soil temperature, thermal properties, ice table
  • Add a routine that compute CO2 + H2O adsorption
  • Minor corrections in PEM.F90

LL

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

Legend:

Unmodified
Added
Removed

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

    r2779 r2794  
    4646  ENDDO
    4747
     48     print *, "jjm_input+1", jjm_input+1
     49     print *, "iim_input+1", iim_input+1
     50     print *, "nslope+1", nslope+1
     51
    4852!  If the difference is too small; there is no evolution
    4953  DO j=1,jjm_input+1
    5054    DO i = 1, iim_input
    5155       DO islope = 1, nslope
     56!         print *, "tendencies_h2o_ice(i,j,islope)LAAA", tendencies_h2o_ice(i,j,islope)
    5257         if(abs(tendencies_h2o_ice(i,j,islope)).LT.1.0E-10) then
    5358            tendencies_h2o_ice(i,j,islope)=0.
    5459         endif
     60!         print *, "tendencies_h2o_ice(i,j,islope)HERE", tendencies_h2o_ice(i,j,islope)
    5561       enddo
    5662    ENDDO
    … …  
    7177  tendencies_h2o_ice_phys(ig0,:) = tendencies_h2o_ice(1,jjm_input+1,:)
    7278
     79!  print *, "tendencies_h2o_ice_physze", tendencies_h2o_ice_phys(:,:)
     80
    7381
    7482END SUBROUTINE compute_tendencies_slope

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

    r2779 r2794  
    1919
    2020!   INPUT
    21   INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points , # subslopes
    22   REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
     21  INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points
     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
    24   REAL,    intent(in) :: ini_surf                      ! Surface of sublimating ice in the GCM output
    25   REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope) ! Boolean to check if it is a sublimating area
    26   REAL,    intent(in) :: global_ave_press_GCM          ! Average pressure in the GCM output
    27   REAL,    intent(in) :: global_ave_press_new          ! Actual pressure
     24  REAL,    intent(in) :: ini_surf
     25  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
     26  REAL,    intent(in) :: global_ave_press_GCM
     27  REAL,    intent(in) :: global_ave_press_new
    2828
    2929
    3030!   OUTPUT
    31   LOGICAL, intent(out) :: STOPPING                     ! Logical : is the criterion reached?
     31  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
    3232
    3333!   local:
    3434!   -----
    35   INTEGER :: i,islope                                  ! Loop
    36   REAL :: present_surf                                 ! Initial/Actual surface of water ice
     35  INTEGER :: i,islope                    ! Loop
     36  REAL :: present_surf  ! Initial/Actual surface of water ice
    3737
    3838!=======================================================================
    … …  
    4141    STOPPING=.FALSE.
    4242
    43 !   computation of the actual surface of sublimating ice
     43!   computation of the actual surface
    4444  present_surf=0.
    4545  do i=1,ngrid
    4646   do islope=1,nslope
    4747      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
    4853         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
    4954      endif
    5055   enddo
    5156  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)
    5264 
    5365!   check of the criterion

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

    r2779 r2794  
    55
    66  USE temps_mod_evol, ONLY: alpha_criterion
    7   use comslope_mod, ONLY: subslope_dist,nslope
     7          use comslope_mod, ONLY: subslope_dist,nslope
    88
    99      IMPLICIT NONE
    … …  
    1919
    2020!   INPUT
    21   INTEGER, intent(in) :: ngrid                         ! # of grid physical grid points
    22   REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
     21  INTEGER, intent(in) :: ngrid                  ! # of grid physical grid points
     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
    24   REAL,    intent(in) :: ini_surf                      ! Surface of sublimating ice in the GCM output
    25   REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope) ! Boolean to check if it is a sublimating area
     24  REAL,    intent(in) :: ini_surf
     25  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
    2626
    2727
    2828!   OUTPUT
    29   LOGICAL, intent(out) :: STOPPING                     ! Logical : is the criterion reached?
     29  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
    3030
    3131!   local:
    3232!   -----
    33   INTEGER :: i,islope                                  ! Loop
    34   REAL :: present_surf                                 ! Initial/Actual surface of water ice
     33  INTEGER :: i,islope                    ! Loop
     34  REAL :: present_surf  ! Initial/Actual surface of water ice
    3535
    3636!=======================================================================
    … …  
    4444    do islope=1, nslope
    4545      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
    4651         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
    4752      endif
    4853    enddo
    4954  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)
    5062 
    5163!   check of the criterion

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

    r2779 r2794  
    2020!   INPUT
    2121
    22   INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope      ! # of grid points along longitude/latitude/ total
    23   REAL, intent(in) ::  cell_area(ngrid)                         ! Area of each cell
     22  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
     24  REAL, intent(in) ::  cell_area(ngrid)
    2425
    2526!   OUTPUT
    26   REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                   ! physical point field : Previous and actual density of water ice
    27   LOGICAL :: STOPPING                                           ! Not used
     27  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                ! physical point field : Previous and actual density of water ice
     28  LOGICAL :: STOPPING
    2829  REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    2930
    … …  
    3334
    3435  INTEGER :: i,j,ig0,islope                                  ! loop variable
     36!  REAL :: not_budget, budget
     37  REAL :: pos_tend, neg_tend, real_coefficient,negative_part
     38  REAL ::  new_tendencies(ngrid)
    3539
    3640  STOPPING=.false.
    3741
    3842
    39 ! Evolution of the co2 ice for each physical point
     43! Evolution of the water ice for each physical point
    4044  do i=1,ngrid
    4145    do islope=1,nslope

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

    r2779 r2794  
    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
    1111
    1212!=======================================================================
    13 ! 
    14 !  Routine that compute the evolution of the water ice.
    15 ! 
    16 !  We suppose that the water exchange uniquely form one point to another
    17 !  (water ice disappear from a sublimating point and increase at another
    18 !  accumulation point). Therfor the total positive and negative tendencies
    19 !  must be equal to conserve the total amount of water
    20 ! 
     13!
     14!  Routine that compute the evolution of the water ice
     15!
    2116!=======================================================================
    2217
    … …  
    2621!   INPUT
    2722
    28   INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope      ! # of grid points along longitude/latitude/ total
    29   REAL, intent(in) ::  cell_area(ngrid)                         ! Area of each cell
     23  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
     25  REAL, intent(in) ::  cell_area(ngrid)
    3026
    3127!   OUTPUT
    32   REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                   ! physical point field : Previous and actual density of water ice
    33   LOGICAL :: STOPPING                                           ! Boolean to alert that there is no sublimating point anymore
     28  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                ! physical point field : Previous and actual density of water ice
     29  LOGICAL :: STOPPING
    3430  REAL, intent(inout) ::  tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    3531
    … …  
    3935
    4036  INTEGER :: i,j,ig0,islope                                  ! loop variable
    41   REAL :: pos_tend, neg_tend, real_coefficient               ! The integral of positive/negative tendencies overall the planet. Ratio
    42   REAL :: negative_part                                      ! After applying the tendencies, some point can be negative, we sum all of them in this variable
    43   REAL :: new_tendencies(ngrid,nslope)                       ! Tendencies after taking into account the non physical negative part above
     37!  REAL :: not_budget, budget
     38  REAL :: pos_tend, neg_tend, real_coefficient,negative_part
     39  REAL ::  new_tendencies(ngrid,nslope)
    4440
    4541
    4642!=======================================================================
    4743
    48 ! Integral of the positive and negative tendencies
     44!  budget=sum(qsurf(:))
     45
    4946  pos_tend=0.
    5047  neg_tend=0.
    … …  
    6259  enddo
    6360
    64 ! These two quantities must be equal to conserve water, therfor we apply a ratio to them
     61!  print *, "pos_tend", pos_tend
     62!  print *, "neg_tend", neg_tend
    6563
    6664  if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then
    … …  
    6866       do islope=1,nslope
    6967       if(tendencies_h2o_ice_phys(i,islope).LT.0) then
     68!          print *, "pos_tend/neg_tend", pos_tend/neg_tend
    7069          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(pos_tend/neg_tend)
    7170       else
    … …  
    7574     enddo
    7675  elseif(neg_tend.LT.pos_tend .and. neg_tend.GT.0) then
     76!          print *, "neg_tend/pos_tend", neg_tend/pos_tend
    7777     do i=1,ngrid
    7878       do islope=1,nslope
    … …  
    8585     enddo
    8686  elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then
    87       call criterion_ice_stop_water_slope(cell_area,1,qsurf(:,:)*0,STOPPING,ngrid,cell_area)
     87
     88!      call criterion_ice_stop(cell_area,1,qsurf*0.,STOPPING,ngrid,cell_area)
     89      call criterion_ice_stop_water_slope(cell_area,1,qsurf(:,:)*0.,STOPPING,ngrid,cell_area)
    8890      do i=1,ngrid
    8991         do islope=1,nslope
    … …  
    9395  endif
    9496
     97 
     98 
     99  negative_part = 0.
     100
    95101
    96102! Evolution of the water ice for each physical point
    97103  do i=1,ngrid
    98104    do islope=1, nslope
     105!      qsurf(i)=qsurf(i)+tendencies_h2o_ice_phys(i)*dt_pem
    99106      qsurf(i,islope)=qsurf(i,islope)+new_tendencies(i,islope)*dt_pem
    100       if (qsurf(i,islope).lt.0) then         ! If we have negative number, we put them to zero add sum them
     107!      budget=budget+tendencies_h2o_ice_phys(i)*dt_pem
     108      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)
    101113        negative_part=negative_part-qsurf(i,islope)*cell_area(i)*subslope_dist(i,islope)
    102114        qsurf(i,islope)=0.
    103115        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
    104124      endif
    105125    enddo
    106126  enddo
    107127
    108 ! This negative part must be put somewhere else to conserve water, therfor we compute new tendencies
     128!  print *, "negative_part", negative_part
    109129  real_coefficient=negative_part/pos_tend
     130!  print *, "real_coefficient", real_coefficient
    110131
    111132  do i=1,ngrid
    … …  
    118139
    119140
     141
     142! Conservation of water ice
     143!  qsurf(:)=qsurf(:)*budget/sum(qsurf(:))
     144
     145
    120146END SUBROUTINE evol_h2o_ice_s_slope

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

    r2779 r2794  
    44! I   Initialisation
    55!    I_a READ run.def , run_pem.def
    6 !    I_b READ starfi_0.nc and start_0.nc
     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
    79
    810! II  Run
    9 !    II_a Compute tendencies
    10 !    II_b Save initial situation
    11 !    II_c Time loop
     11!    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
    1215
    1316! III Output
    14 !    III_a Write startfi.nc
     17!    III_a Recomp tendencies for the start and startfi
     18!    III_b Write start and starfi.nc
     19!    III_c Write start_pem
    1520
    1621!------------------------
    … …  
    2126!module needed for INITIALISATION
    2227      use phyetat0_mod, only: phyetat0
    23       use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer
     28      use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
    2429      use surfdat_h, only: tsurf, co2ice, emis,&
    2530      &                    qsurf,watercap, ini_surfdat_h, &
    2631                           albedodat, zmea, zstd, zsig, zgam, zthe, &
    27       &                    hmons, summit, base
     32                           hmons, summit, base,albedo_h2o_frost, &
     33                          frost_albedo_threshold,emissiv
    2834      use dimradmars_mod, only: totcloudfrac, albedo, &
    2935                                ini_dimradmars_mod
    3036      use turb_mod, only: q2, wstar, ini_turb_mod
    3137      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
    3240      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
    3341                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
    3442                        nf90_inquire_dimension,nf90_close
    3543      use phyredem, only: physdem0, physdem1
    36       use tracer_mod, only: noms,nqmx ! tracer names
     44      use tracer_mod, only: noms,nqmx,igcm_h2o_ice ! tracer names
    3745
    3846! For phyredem :
    … …  
    4957      USE infotrac
    5058      USE temps_mod_evol, ONLY: nyear, dt_pem
    51       USE vertical_layers_mod, ONLY: ap,bp
    52 
    53       use comslope_mod, ONLY: nslope,def_slope,def_slope_mean, &
     59!     USE vertical_layers_mod, ONLY: ap,bp
     60      USE comcstfi_h, only: pi
     61
     62      USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, &
    5463                           subslope_dist,co2ice_slope, &
    5564                           tsurf_slope,tsoil_slope,fluxgrd_slope,&
    … …  
    5968                           iflat
    6069
    61 
     70      USE geometry_mod, only: longitude_deg,latitude_deg
     71
     72      use pemredem, only:  pemdem1
     73      USE soil_evolution_mod, ONLY: soil_pem_CN
     74!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
     75
     76      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
     82      use adsorption_mod, only : regolith_co2adsorption
    6283  IMPLICIT NONE
    6384
    6485  include "dimensions.h"
    6586  include "paramet.h"
     87
     88  INTEGER ngridmx
     89  PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
     90
    6691  include "comdissnew.h"
    6792  include "comgeom.h"
    6893  include "iniprint.h"
    69 
    7094! Same variable's name as in the GCM
    7195
    72       INTEGER :: ngrid      !Number of physical grid points
    73       INTEGER :: nlayer     !Number of vertical layer
    74       INTEGER :: nq         !Number of tracer
    75       INTEGER :: day_ini    !First day of the simulation
    76       REAL :: pday          !Physical day
    77       REAL :: time_phys     !Same as GCM
    78       REAL :: ptimestep     !Same as GCM
    79       REAL :: ztime_fin     !Same as GCM
     96  INTEGER :: ngrid      !Number of physical grid points
     97  INTEGER :: nlayer     !Number of vertical layer
     98  INTEGER :: nq         !Number of tracer
     99  INTEGER :: day_ini    !First day of the simulation
     100  REAL :: pday          !Physical day
     101  REAL :: time_phys     !Same as GCM
     102  REAL :: ptimestep     !Same as GCM
     103  REAL :: ztime_fin     !Same as GCM
    80104
    81105! Variable for reading start.nc
    82       character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM
     106      character (len = *), parameter :: FILE_NAME_start = "start_0.nc" !Name of the file used for initialsing the PEM
    83107  !   variables dynamiques
    84108  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    … …  
    86110  REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
    87111  REAL ps(ip1jmp1)                       ! pression  au sol
    88   REAL, dimension(:),allocatable :: ps_phys !(ngrid) pression  au sol
     112  REAL, dimension(:),allocatable :: ps_phys !(ngrid)                       ! pression  au sol
     113  REAL, dimension(:,:),allocatable :: ps_phys_timeseries !(ngrid x timelen) ! pression  au sol instantannées
     114  REAL, dimension(:,:),allocatable :: ps_phys_timeseries_yr1 !(ngrid x timelen) ! pression  au sol instantannées for the first year of the gcm
     115
    89116  REAL masse(ip1jmp1,llm)                ! masse d'air
    90117  REAL phis(ip1jmp1)                     ! geopotentiel au sol
    … …  
    93120! Variable for reading starfi.nc
    94121
    95       character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
     122      character (len = *), parameter :: FILE_NAME = "startfi_0.nc" !Name of the file used for initialsing the PEM
    96123      integer :: ncid, varid,status                                !Variable for handling opening of files
    97124      integer :: phydimid, subdimid, nlayerdimid, nqdimid          !Variable ID for Netcdf files
    98125      integer :: lonvarid, latvarid, areavarid,sdvarid             !Variable ID for Netcdf files
     126      integer :: apvarid,bpvarid                                   !Variable ID for Netcdf files
    99127
    100128! Variable for reading starfi.nc and writting restartfi.nc
    … …  
    102130      REAL, dimension(:),allocatable :: longitude                  !Longitude read in FILE_NAME and written in restartfi
    103131      REAL, dimension(:),allocatable  :: latitude                  !Latitude read in FILE_NAME and written in restartfi
     132      REAL, dimension(:),allocatable :: ap                         !Coefficient ap read in FILE_NAME_start and written in restart
     133      REAL, dimension(:),allocatable :: bp                         !Coefficient bp read in FILE_NAME_start and written in restart
    104134      REAL, dimension(:),allocatable  :: cell_area                 !Cell_area read in FILE_NAME and written in restartfi
    105135      REAL :: Total_surface                                        !Total surface of the planet
    … …  
    113143      REAL, dimension(:),allocatable  :: tendencies_co2_ice_phys   ! physical point field : Tendency of evolution of perenial co2 ice over a year
    114144
    115 ! Variable for initial situation of water ice
    116 
    117145      REAL :: ini_surf                                             ! Initial surface of sublimating water ice
    118       REAL :: ini_surf_h2o                                         ! Initial surface of sublimating water ice
     146      REAL :: ini_surf_h2o                                             ! Initial surface of sublimating water ice
    119147      REAL, dimension(:),allocatable  :: initial_h2o_ice           ! physical point field : Logical array indicating sublimating point
    120 
    121 ! Variable for initial situation of co2 ice
    122148
    123149      REAL :: ini_surf_co2                                         ! Initial surface of sublimating co2 ice
    124150      REAL, dimension(:),allocatable  :: initial_co2_ice           ! physical point field : Logical array indicating sublimating point of co2 ice
    125 
    126 ! Variable for needed to compute tendencies for water ice
    127151
    128152      REAL , dimension(:,:), allocatable :: min_h2o_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
    … …  
    130154      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
    131155
    132 ! Variable for needed to compute tendencies for co2 ice
    133 
    134156      REAL , dimension(:,:), allocatable :: min_co2_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
    135157      REAL , dimension(:,:), allocatable :: min_co2_ice_s_2        ! LON x LAT field : minimum of water ice at each point for the second year
    136158      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
    137159
    138       REAL :: global_ave_press_GCM           ! physical point field : Global average pressure in the GCM output
     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
     162
     163      REAL :: global_ave_press_GCM
    139164      REAL :: global_ave_press_old           ! physical point field : Global average pressure of initial/previous time step
    140165      REAL :: global_ave_press_new           ! physical point field : Global average pressure of current time step
    141166
    142       REAL , dimension(:,:), allocatable :: ZPLEV_new ! ngrid x nlayer+1 : Pressure level at the end of the iteration
    143       REAL , dimension(:,:), allocatable :: ZPLEV_OLD ! ngrid x nlayer+1 : Pressure level at the beginning of the iteration
    144 
    145       INTEGER :: year_iter        !Counter for the number of PEM iteration
    146       INTEGER :: year             ! Year read in the startfi files and written back in it counting the year iterations
    147 
     167      REAL , dimension(:,:), allocatable ::  zplev_new
     168      REAL , dimension(:,:), allocatable :: zplev_old
     169      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
    148179      LOGICAL :: STOPPING_water   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
    149180      LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
    150 
    151       LOGICAL :: STOPPING_co2     ! Logical : is the criterion (% of change in the surface of sublimating co2 ice) reached?
    152       LOGICAL :: STOPPING_1_co2   ! Logical : is there still co2 ice to sublimate?
    153 
    154       REAL, dimension(:),allocatable  :: q_co2_GCM ! Initial amount of co2 in the first layer
    155       REAL ps_GCM(ip1jmp1)                         ! pression  au sol donné par le GCM
    156 
    157       real ,allocatable :: vmr_co2_gcm_phys(:,:)   !(ngrid x Timelen) : co2 volume mixing ratio given by the GCM
    158       REAL, ALLOCATABLE ::  vmr_co2_gcm(:,:,:)     ! iim+1 x jjm+1 x Timelen : co2 volume mixing ratio given by the GCM
    159 
    160       REAL, ALLOCATABLE ::  ps_GCM_1(:,:,:)        ! iim+1 x jjm+1 x Timelen : Surface pressure in the first GCM year
    161       REAL, ALLOCATABLE ::  ps_GCM_2(:,:,:)        ! iim+1 x jjm+1 x Timelen : Surface pressure in the second GCM year
    162 
    163       integer :: timelen                           ! Number of time point in the GCM diagfi output
    164 
    165       REAL, ALLOCATABLE :: p(:,:)                  !(ngrid,llmp1) : Pressure
     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?
     183
     184      REAL, dimension(:,:,:),allocatable  :: q_co2_GCM ! Initial amount of co2 in the first layer
     185
     186      real,save :: m_co2, m_noco2, A , B, mmean
     187      real ,allocatable :: vmr_co2_gcm_phys(:,:) !(ngrid) ! co2 volume mixing ratio
     188      real ,allocatable :: vmr_co2_pem_phys(:,:) !(ngrid) ! co2 volume mixing ratio
     189      real ,allocatable :: q_h2o_GCM_phys(:,:)
     190      real ,allocatable :: q_co2_GCM_phys(:,:)
     191      real ,allocatable :: q_co2_PEM_phys(:,:)
     192      REAL, ALLOCATABLE ::  ps_GCM(:,:,:)
     193      REAL, ALLOCATABLE :: ps_GCM_yr1(:,:,:)
     194      REAL, ALLOCATABLE ::  vmr_co2_gcm(:,:,:)
     195      REAL, ALLOCATABLE ::  q_h2o_GCM(:,:,:)
     196      REAL ,allocatable ::  q_h2o_PEM_phys(:,:)
     197      real ,allocatable :: q_phys(:,:,:)
     198      integer :: timelen
     199      REAL :: ave
     200
     201      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
    166208
    167209!!!!!!!!!!!!!!!!!!!!!!!! SLOPE
    168       REAL ,allocatable :: watercap_slope(:,:)                           ! (ngrid,nslope)
    169       REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of co2 ice at each point for the first year
    170       REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of co2 ice at each point for the second year
     210      character*2 :: str2
     211      REAL ,allocatable :: watercap_slope(:,:)    !(ngrid,nslope)
     212      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of water ice at each point for the first year
     213      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of water ice at each point for the second year
    171214      REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_1        ! LON x LAT field : minimum of water ice at each point for the first year
    172215      REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_2        ! LON x LAT field : minimum of water ice at each point for the second year
    173216
     217
     218      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
     223      REAL, dimension(:,:),allocatable  :: initial_h2o_ice_slope           ! physical point field : Logical array indicating sublimating point of h2o ice
    174224      REAL, dimension(:,:),allocatable  :: initial_co2_ice_slope           ! physical point field : Logical array indicating sublimating point of co2 ice
    175       REAL, dimension(:,:),allocatable  :: initial_h2o_ice_slope           ! physical point field : Logical array indicating sublimating point of co2 ice
    176225
    177226      REAL , dimension(:,:,:), allocatable :: tendencies_co2_ice_slope     ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
    178227      REAL , dimension(:,:,:), allocatable :: tendencies_h2o_ice_slope     ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
    179228      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope   ! physical point field : Tendency of evolution of perenial co2 ice over a year
     229      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope_ini ! physical point field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
    180230      REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical point field : Tendency of evolution of perenial co2 ice over a year
    181231
    182232      REAL, PARAMETER :: co2_hmax = 10                    ! Maximum height  for CO2 deposit on slopes (m)
    183       REAL, PARAMETER :: rho_co2 = 1600                   ! CO2 ice density (kg/m^3)
    184       INTEGER :: iaval                                    ! Index of the neighboord slope ()
     233      REAL, PARAMETER :: rho_co2 = 1600           ! CO2 ice density (kg/m^3)
     234      INTEGER :: iaval                            ! Index of the neighboord slope ()
    185235      REAL , dimension(:,:), allocatable :: flag_co2flow(:,:)   !(ngrid,nslope)          ! To flag where there is a glacier flow
    186236      REAL , dimension(:), allocatable :: flag_co2flow_mesh(:)  !(ngrid)          ! To flag where there is a glacier flow
    187237
    188238
    189 
     239!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL
     240
     241
     242     REAL, ALLOCATABLE :: tsurf_ave(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
     243     REAL, ALLOCATABLE  :: tsurf_ave_phys(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature [K]
     244     REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
     245     REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
     246     REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) !IG x SLOPE field : Averaged Soil Temperature [K]
     247
     248     REAL, ALLOCATABLE :: TI_GCM(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Thermal Inertia  [SI]
     249     REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:) !LONX LAT x SLOPE XTULES field : NOn averaged Surf Temperature [K]
     250     REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) !IG x SLOPE XTULES field : NOn averaged Surf Temperature [K]
     251
     252     REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
     253     REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
     254
     255     REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
     256     REAL, ALLOCATABLE  :: tsurf_ave_phys_yr1(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
     257
     258
     259
     260     REAL :: tol_soil = 1.e-6
     261     INTEGER :: countmax = 15
     262     INTEGER :: countloop
     263     REAL, ALLOCATABLE :: tsoil_phys_PEM_saved(:,:,:)
     264     REAL :: error_tsoil
     265
     266
     267     LOGICAL :: firstcall
     268!    REAL, PARAMETER :: daysec=88775.              ! duree du sol (s)  ~88775 s
     269     REAL, PARAMETER :: year_day = 669
     270     REAL, PARAMETER :: year_step = 1
     271     REAL :: timestep
     272
     273     REAL, ALLOCATABLE :: inertiesoil(:,:) ! Thermal inertia of the mesh for restart
     274
     275     REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:) ! Averaged GCM Thermal Inertia  [SI]
     276     REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Averaged GCM Thermal Inertia  [SI]
     277
     278     REAL,ALLOCATABLE  :: q_h2o_PEM_phys_ave(:) ! averaged water vapor content
     279     REAL,ALLOCATABLE  :: interp_coef(:)
     280
     281     REAL,ALLOCATABLE  :: ice_depth(:,:)
     282     REAL,ALLOCATABLE  :: TI_locslope(:,:)
     283     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:)
     284     REAL,ALLOCATABLE  :: Tsoil_locslope_CN(:,:)
     285     REAL,ALLOCATABLE  :: Tsurf_locslope(:)
     286     REAL,ALLOCATABLE  :: Tsurf_locslope_prev(:)
     287     REAL,ALLOCATABLE  :: alph_locslope(:,:)
     288     REAL,ALLOCATABLE  :: beta_locslope(:,:)   
     289     REAL :: kcond
     290     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)
     291     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)
     292     REAL :: totmass_adsorbded
    190293!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    191294
    192295! Loop variable
    193      INTEGER :: i,j,ig0,l,ig,nnq,t,islope
    194 ! Constant
    195      REAL :: pi,beta,alpha
    196 
     296     LOGICAL :: bool_sublim
     297     INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop
     298     REAL :: beta = 3182.48
     299     REAL :: alpha = 23.3494
     300 
    197301! Parallel variables
    198302      is_sequential=.true.
    … …  
    205309      time_phys=0. !test
    206310
    207 !------------------------
    208 
    209 !=======================================================
     311!      n_band_lat=18
     312
     313      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
     314      m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
     315      A =(1/m_co2 - 1/m_noco2)
     316      B=1/m_noco2
     317
     318
     319
     320!------------------------
    210321
    211322! I   Initialisation
    … …  
    219330      CALL conf_evol( 99, .TRUE. )
    220331
     332
    221333!------------------------
    222334
    … …  
    262374
    263375!----------------------------READ start.nc ---------------------
    264 
    265     call infotrac_init
     376   
     377     call infotrac_init
    266378
    267379     allocate(q(ip1jmp1,llm,nqtot))
    268 
    269380     CALL dynetat0(FILE_NAME_start,vcov,ucov, &
    270381                    teta,q,masse,ps,phis, time_0)
    271382
    272      CALL iniconst
     383     CALL iniconst !new
    273384     CALL inigeom
     385     allocate(ap(nlayer+1))
     386     allocate(bp(nlayer+1))
     387     status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid)
     388     status = nf90_inq_varid(ncid, "ap", apvarid)
     389     status = nf90_get_var(ncid, apvarid, ap)
     390     status = nf90_inq_varid(ncid, "bp", bpvarid)
     391     status = nf90_get_var(ncid, bpvarid, bp)
     392     status =nf90_close(ncid)
     393     
     394       
     395   
     396    daysec=88775.
     397    dtphys=0
     398
     399 
    274400
    275401    CALL iniphysiq(iim,jjm,llm, &
    … …  
    278404          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
    279405          iflag_phys)
    280 
    281      year=day_ini
    282 
    283 
     406   DO nnq=1,nqtot
     407    if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
     408   ENDDO
     409
     410 
    284411!----------------------------reading ---------------------
    285412
    … …  
    287414! First we read the initial state (starfi.nc)
    288415
     416
    289417       allocate(watercap_slope(ngrid,nslope))
     418       allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
     419       allocate(q_h2o_PEM_phys_ave(ngrid))
     420      allocate(inertiesoil(ngrid,nsoilmx))
     421
    290422
    291423         CALL phyetat0 (FILE_NAME,0,0, &
    … …  
    294426              tsurf,tsoil,albedo,emis,   &
    295427              q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar,     &
    296               watercap,nslope,tsurf_slope,                       &
     428              watercap,inertiesoil,nslope,tsurf_slope,           &
    297429              tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
    298               subslope_dist,albedo_slope,emiss_slope,            &
     430              subslope_dist,albedo_slope,emiss_slope, TI_GCM_start,     &
    299431              qsurf_slope,watercap_slope)
    300432
    301 ! Slope : Definition of the flat slope
     433
     434
     435     deallocate(TI_GCM_start) !not used then
     436
     437     DO i=1,ngrid
     438         DO nnq=1,nqtot
     439           DO islope=1,nslope
     440             if(qsurf_slope(i,nnq,islope).LT.0) then
     441                qsurf_slope(i,nnq,islope)=0.
     442             endif
     443           enddo
     444         enddo
     445       enddo
     446
     447!     Define some slope statistics
     448
     449
     450     
    302451
    303452       iflat=1
    … …  
    306455           abs(def_slope_mean(iflat)))THEN
    307456           iflat = islope
    308          ENDIF
     457         ENDIF     
     458
    309459       ENDDO
    310460       PRINT*,'Flat slope for islope = ',iflat
    311461       PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
    312462
    313      allocate(flag_co2flow(ngrid,nslope))
    314      allocate(flag_co2flow_mesh(ngrid))
    315 
    316463       flag_co2flow(:,:) = 0.     
    317464       flag_co2flow_mesh(:) = 0.
    318465
    319 
    320 ! Then we read the evolution of water and co2 ice over the first year of the GCM run, saving only the minimum value
    321 
     466!------------------------
     467
     468! I   Initialisation
     469!    I_a READ run.def , run_pem.def
     470!    I_b READ starfi_0.nc
     471!    I_c READ GCM data and convert to the physical grid
     472
     473!------------------------
     474
     475! 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
     476
     477
     478     call nb_time_step_GCM("concat_year_one.nc",timelen)
    322479     allocate(min_h2o_ice_s_1(iim+1,jjm+1))
    323480     allocate(min_co2_ice_s_1(iim+1,jjm+1))
    324 
    325      allocate(ps_GCM_1(iim+1,jjm+1,2676))
    326 
     481     allocate(vmr_co2_gcm(iim+1,jjm+1,timelen))
     482     allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
     483     allocate(q_co2_GCM(iim+1,jjm+1,timelen))
     484     allocate(ps_GCM(iim+1,jjm+1,timelen))
     485     allocate(ps_GCM_yr1(iim+1,jjm+1,timelen))
    327486     allocate(min_co2_ice_slope_1(iim+1,jjm+1,nslope))
    328487     allocate(min_h2o_ice_slope_1(iim+1,jjm+1,nslope))
    329 
    330      allocate(vmr_co2_gcm(iim+1,jjm+1,2676))
    331 
    332      call pemetat0 ("concat_year_one.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_1,timelen,min_co2_ice_slope_1,min_h2o_ice_slope_1,nslope)
    333 
    334 ! Then we read the evolution of water and co2 ice over the second year of the GCM run, saving only the minimum value
     488     allocate(tsurf_ave(iim+1,jjm+1,nslope))
     489     allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope))
     490     allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
     491     allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
     492     allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
     493     allocate(tsurf_ave_phys(ngrid,nslope))
     494     allocate(tsurf_ave_phys_yr1(ngrid,nslope))
     495     allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
     496     allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
     497     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
     498
     499     allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen))
     500     allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen))
     501     allocate(co2_ice_GCM_phys_slope(ngrid,nslope,timelen))
     502     allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen))
     503     allocate(Tsurfave_before_saved(ngrid,nslope))
     504     allocate(tsoil_phys_PEM_saved(ngrid,nsoilmx_PEM,nslope))
     505     allocate(delta_co2_adsorbded(ngrid))
     506     allocate(q_phys(ngrid,llm,nqtot))
     507
     508
     509
     510     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,&   
     511                       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)
     512
     513
     514! 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
    335515
    336516     allocate(min_h2o_ice_s_2(iim+1,jjm+1))
    337517     allocate(min_co2_ice_s_2(iim+1,jjm+1))
    338 
    339518     allocate(min_co2_ice_slope_2(iim+1,jjm+1,nslope))
    340519     allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
    341520
    342      allocate(ps_GCM_2(iim+1,jjm+1,2676))
    343 
    344      call pemetat0 ("concat_year_two.nc", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_2,timelen,min_co2_ice_slope_2,min_h2o_ice_slope_2,nslope)
    345 
    346 ! We read the evolutaion of water ice over the first year of the GCM run, saving the first and last state
     521     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, &
     522                  nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope)
     523
     524
     525
     526
     527! The variables in the dynamic grid are transfered to the physical grid
     528
    347529
    348530     allocate(vmr_co2_gcm_phys(ngrid,timelen))
    349 
    350   vmr_co2_gcm_phys(1,:)=vmr_co2_gcm(1,1,:)
    351 
    352   ig0 = 2
    353   DO j=2,jjm
    354     DO i = 1, iim
    355        vmr_co2_gcm_phys(ig0,:)  =vmr_co2_gcm(i,j,:)
    356        ig0= ig0 + 1
    357     ENDDO
    358   ENDDO
    359 
    360   vmr_co2_gcm_phys(ig0,:) = vmr_co2_gcm(1,jjm+1,:)
    361 
    362 !=======================================================
    363 
    364 ! II  Run
    365 !    II_a Compute tendencies
    366 
    367 !------------------------
    368 
    369 !---------------------------- RUN ---------------------
     531     allocate(vmr_co2_pem_phys(ngrid,timelen))
     532     allocate(q_h2o_GCM_phys(ngrid,timelen))
     533     allocate(q_h2o_PEM_phys(ngrid,timelen))
     534     allocate(q_co2_GCM_phys(ngrid,timelen))
     535     allocate(q_co2_PEM_phys(ngrid,timelen))
     536
     537
     538     CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys)
     539     CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys)
     540     CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys)
     541
     542     do nnq=1,nqtot
     543
     544        call gr_dyn_fi(llm,iip1,jjp1,ngridmx,q(:,:,nnq),q_phys(:,:,nnq))
     545     enddo
     546
     547
     548     deallocate(vmr_co2_gcm)
     549     deallocate(q_h2o_GCM)
     550     deallocate(q_co2_GCM)
     551
     552
     553    q_co2_PEM_phys(:,:)=  q_co2_GCM_phys(:,:)
     554    q_h2o_PEM_phys(:,:)=  q_h2o_GCM_phys(:,:)
     555
     556
     557     allocate(ps_phys(ngrid))
     558     allocate(ps_phys_timeseries(ngrid,timelen))
     559     allocate(ps_phys_timeseries_yr1(ngrid,timelen))
     560
     561
     562
     563
     564
     565     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
     566
     567
     568! for the time series:
     569
     570
     571
     572    call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
     573    call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
     574    deallocate(ps_GCM)
     575    deallocate(ps_GCM_yr1)
     576
     577
     578!!!!!!!!!!!!!!!!!!!!! Initialisation for soil_PEM
     579
     580      call end_comsoil_h_PEM
     581      call ini_comsoil_h_PEM(ngrid,nslope)
     582      timestep=year_day*daysec/year_step
     583
     584
     585      allocate(ice_depth(ngrid,nslope))
     586      allocate(TI_GCM_phys(ngrid,nsoilmx,nslope))
     587
     588      DO islope = 1,nslope
     589          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope))
     590      ENDDO
     591
     592
     593      call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
     594
     595      deallocate(TI_GCM)
     596
     597
     598      DO islope = 1,nslope
     599
     600        DO l=1,nsoilmx
     601           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
     602        ENDDO
     603        DO l=nsoilmx+1,nsoilmx_PEM
     604           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
     605        ENDDO
     606      ENDDO
     607
     608
     609
     610       deallocate(tsoil_ave_yr1)
     611
     612
     613
     614
     615
     616
     617
     618
     619
     620
     621      DO islope = 1,nslope
     622
     623        DO l=1,nsoilmx
     624           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
     625        ENDDO
     626        DO l=nsoilmx+1,nsoilmx_PEM
     627           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
     628        ENDDO
     629      ENDDO
     630     
     631
     632
     633       deallocate(tsoil_ave)
     634
     635
     636      DO islope = 1,nslope
     637       
     638        DO l=1,nsoilmx
     639          DO t=1,timelen
     640           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
     641          ENDDO
     642        ENDDO
     643      ENDDO
     644
     645      deallocate(tsoil_GCM_timeseries)
     646     
     647 
     648
     649
     650       CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
     651       
     652       CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
     653
     654       deallocate(tsurf_ave)
     655       deallocate(tsurf_ave_yr1)
     656
     657   
     658
     659
     660      DO islope = 1,nslope
     661
     662        DO t=1,timelen
     663          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
     664        enddo
     665      enddo
     666      deallocate(tsurf_GCM_timeseries)
     667
     668
     669  DO islope = 1,nslope
     670
     671        DO t=1,timelen
     672          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
     673        enddo
     674      enddo
     675
     676
     677
     678      deallocate(co2_ice_GCM_slope)
     679
     680
     681
     682!------------------------
     683
     684! I   Initialisation
     685!    I_a READ run.def , run_pem.def
     686!    I_b READ starfi_0.nc
     687!    I_c READ GCM data and convert to the physical grid
     688!    I_d Compute tendencies & Save initial situation
     689
     690!------------------------
     691
     692
     693
     694
    370695
    371696!----- Compute tendencies
    … …  
    380705     allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
    381706     allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
     707     allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope))
     708
    382709
    383710     allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope))
    … …  
    389716             min_h2o_ice_s_2,iim,jjm,ngrid,tendencies_h2o_ice_phys)
    390717
    391 !  Compute the tendencies of the evolution of co2 ice over the years
    392 
    393718      call compute_tendencies(tendencies_co2_ice,min_co2_ice_s_1,&
    394719             min_co2_ice_s_2,iim,jjm,ngrid,tendencies_co2_ice_phys)
    395720
    396 !  Compute the tendencies of the evolution of water ice over the years with slopes
    397 
    398721      call compute_tendencies_slope(tendencies_co2_ice_slope,min_co2_ice_slope_1,&
    399722             min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope)
    400723
    401 !  Compute the tendencies of the evolution of co2 ice over the years with slopes
     724
     725      tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:)
    402726
    403727      call compute_tendencies_slope(tendencies_h2o_ice_slope,min_h2o_ice_slope_1,&
    404728             min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope)
    405729
    406 !------------------------
    407 
    408 ! II  Run
    409 !    II_a Compute tendencies
    410 !    II_b Save initial situation
    411 
    412 !------------------------
    413 
     730     
    414731!----- Save initial situation
    415732
    416733     allocate(initial_h2o_ice(ngrid))
    417 
    418734     allocate(initial_co2_ice(ngrid))
    419      allocate(q_co2_GCM(ngrid))
    420 
     735     allocate(initial_co2_ice_sublim_slope(ngrid,nslope))
    421736     allocate(initial_co2_ice_slope(ngrid,nslope))
    422737     allocate(initial_h2o_ice_slope(ngrid,nslope))
    423 
    424 ! We save the places where water/co2 ice is sublimating
     738     year_iter=0
     739
     740! We save the places where water ice is sublimating
    425741  do i=1,ngrid
    426742      if (tendencies_h2o_ice_phys(i).LT.0) then
    … …  
    432748
    433749      if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
     750         initial_co2_ice_sublim_slope(i,islope)=1.
     751      else
     752         initial_co2_ice_sublim_slope(i,islope)=0.         
     753      endif
     754      if (co2ice_slope(i,islope).GT.0) then
    434755         initial_co2_ice_slope(i,islope)=1.
    435756      else
    … …  
    449770  ini_surf_h2o=0.
    450771  Total_surface=0.
    451 
    452772  do i=1,ngrid
    453773    if (initial_h2o_ice(i).GT.0.5) then
    … …  
    455775    endif
    456776    do islope=1,nslope
    457       if (initial_co2_ice_slope(i,islope).GT.0.5) then
     777      if (initial_co2_ice_sublim_slope(i,islope).GT.0.5) then
    458778         ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
    459779      endif
    … …  
    465785  enddo
    466786
    467 ! We put the pressure on the physical grid
    468 
    469   allocate(ps_phys(ngrid))
    470 
    471   ig0 = iim+1
    472   ps_phys(1)=ps(ig0)
    473   ig0=ig0+1
    474 
    475   DO i = 2, ngrid-1
    476      if(modulo(ig0,iim+1).eq.0 .and. i.gt.3) then
    477        ig0=ig0+1
    478      endif
    479      ps_phys(i) = ps(ig0)
    480      ig0= ig0 + 1
    481   ENDDO
    482 
    483   ps_phys(ngrid)=ps(ig0)
    484 
    485 ! We compute the global average pressure of the initial state
    486 
    487   global_ave_press_old=0.
    488   do i=1,ngrid
    489     global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface
    490   enddo
    491 
    492   global_ave_press_GCM=global_ave_press_old
    493 
    494 ! We compute the initial pressure level
    495 
     787     print *, "ini_surf_co2=", ini_surf_co2
     788     print *, "ini_surf=", ini_surf
     789     print *, "ini_surf_h2o=", ini_surf_h2o
     790
     791
     792     allocate(local_old_press(ngrid))
     793     allocate(local_new_press(ngrid))
     794   
    496795     allocate(zplev_new(ngrid,nlayer+1))
    497796     allocate(zplev_old(ngrid,nlayer+1))
    498797
     798
     799     global_ave_press_old=0.
     800     do i=1,ngrid
     801       global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface
     802     enddo
     803
     804     global_ave_press_GCM=global_ave_press_old
     805     print *, "global_ave_press_old", global_ave_press_old
     806     
     807
     808
     809!----- Time loop
     810     print *, "Before timeloop"
     811     allocate(flag_co2flow(ngrid,nslope))
     812     allocate(flag_co2flow_mesh(ngrid))
     813
     814     firstcall=.TRUE.
     815
     816
     817
     818
     819
     820!------------------------
     821
     822! I   Initialisation
     823!    I_a READ run.def , run_pem.def
     824!    I_b READ starfi_0.nc
     825!    I_c READ GCM data and convert to the physical grid
     826!    I_d Compute tendencies & Save initial situation
     827!    I_e Read startfi_PEM
     828
     829!------------------------
     830
     831
     832! now we complete with the PEMstart
     833
     834      call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
     835      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,&
     836      tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), co2_adsorbded_phys,delta_co2_adsorbded)
     837totmass_adsorbded = 0.
     838     do ig = 1,ngrid
     839      do islope =1, nslope
     840       do l = 1,nsoilmx_PEM - 1
     841
     842          totmass_adsorbded = totmass_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     843          subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
     844          cell_area(ig)
     845       enddo
     846
     847      enddo
     848
     849     enddo
     850
     851     write(*,*) "tot mass in the regolith", totmass_adsorbded
     852     stop
     853     deallocate(tsurf_ave_phys_yr1)
     854     deallocate(ps_phys_timeseries_yr1)
     855     deallocate(tsoil_ave_phys_yr1)   
     856
     857
     858!--------------------------- END INITIALISATION ---------------------
     859
     860
     861
     862
     863!---------------------------- RUN ---------------------
     864
     865
     866
     867
     868!------------------------
     869
     870! II  Run
     871!    II_a update pressure,ice and tracers
     872
     873!------------------------
     874
     875   
     876     do while (.true.)
     877! II.a.1. global pressure
     878     global_ave_press_new=global_ave_press_old
     879
     880     do i=1,ngrid
     881       do islope=1,nslope
     882           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
     883       enddo
     884     
     885      global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
     886     enddo
     887! II.a.2. tracers
     888     allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
    499889        DO l=1,nlayer+1
    500          DO ig=1,ip1jmp1
    501           zplev_old(ig,l) = ap(l)  + bp(l)*ps(ig)
     890         DO ig=1,ngrid
     891          zplev_old(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
     892          zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
    502893         ENDDO
    503894        ENDDO
    504895
    505 
    506 !=======================================================
     896     do i=1,ip1jmp1
     897
     898       ps(i)=ps(i)*global_ave_press_new/global_ave_press_old
     899
     900     enddo
     901     do i = 1,ngrid
     902       ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_old
     903       ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old
     904     enddo
     905
     906     allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
     907
     908        do l=1,nlayer+1
     909         do ig=1,ngrid
     910          zplev_new(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
     911          zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
     912         enddo
     913        enddo
     914
     915
     916      print *, "1 is okay"
     917
     918        DO nnq=1,nqtot
     919        if (noms(nnq).NE."co2") then
     920          DO l=1,llm-1
     921            DO ig=1,ngrid
     922               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))
     923            ENDDO
     924            q(:,llm,nnq)=q(:,llm-1,nnq)
     925          ENDDO
     926        else
     927          DO l=1,llm-1
     928            DO ig=1,ngrid
     929                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))  &
     930                                + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
     931                                      (zplev_old(ig,l)-zplev_old(ig,l+1))     )  / &
     932                                      (zplev_new(ig,l)-zplev_new(ig,l+1))
     933            ENDDO
     934           q(:,llm,nnq)=q(:,llm-1,nnq)
     935          ENDDO
     936        endif
     937        ENDDO
     938
     939
     940 DO nnq=1,nqtot
     941       DO ig=1,ngrid
     942         DO l=1,llm-1
     943          if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number" .or. noms(nnq).NE."ccn_number") ) then
     944              extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
     945              q(ig,l,nnq)=1.
     946              q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
     947          endif
     948          if(q(ig,l,nnq).LT.0) then
     949              q(ig,l,nnq)=1E-30
     950          endif
     951         ENDDO
     952       enddo
     953   enddo
     954
     955
     956
     957
     958  l=1
     959  DO ig=1,ngrid
     960      DO t = 1, timelen
     961
     962         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))
     963         if(q_h2o_PEM_phys(ig,t).LT.0) then
     964           q_h2o_PEM_phys(ig,t)=1E-30
     965         endif
     966         if(q_h2o_PEM_phys(ig,t).GT.1) then
     967           q_h2o_PEM_phys(ig,t)=1.
     968         endif
     969
     970   enddo
     971  enddo
     972
     973
     974
     975  DO ig=1,ngrid
     976      DO t = 1, timelen
     977         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))  &
     978                                + (   (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))  -       &
     979                                      (zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))     )  / &
     980                                      (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
     981         if (q_co2_PEM_phys(ig,t).LT.0) then
     982              q_co2_PEM_phys(ig,t)=1E-30
     983         elseif (q_co2_PEM_phys(ig,t).GT.1) then
     984              q_co2_PEM_phys(ig,t)=1.
     985         endif
     986         mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
     987         vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
     988      ENDDO
     989  ENDDO
     990   
     991     deallocate(zplev_new_timeseries)
     992     deallocate(zplev_old_timeseries)
     993
     994
     995! II.a.3. Evolution of the ice
     996
     997
     998     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)
     999      print *, "4 is okay"
     1000
     1001     call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
     1002      print *, "5 is okay"
     1003
    5071004
    5081005!------------------------
    5091006
    5101007! II  Run
    511 !    II_a Compute tendencies
    512 !    II_b Save initial situation
    513 !    II_c Time loop
    514 
    515 !------------------------
    516 
    517 
    518 !----- Time loop
    519 
    520      year_iter=0
    521 
    522      do while (.true.)
    523 
    524 ! Compute the new global average pressure given the tendencies of co2 condensation
    525      global_ave_press_new=global_ave_press_old
    526      do i=1,ngrid
    527        do islope=1,nslope
    528            global_ave_press_new=global_ave_press_new-dt_pem*g*cell_area(i)*tendencies_co2_ice_phys_slope(i,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
    529        enddo
    530      enddo
    531 
    532 ! Compute the new surface pressure at each physical point
    533      do i=1,ngrid
    534        ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_old
    535      enddo
    536 
    537 ! Compute the new surface pressure at each dynamical point
    538      do i=1,ip1jmp1
    539        ps(i)=ps(i)*global_ave_press_new/global_ave_press_old
    540      enddo
    541 
    542 ! Compute the new pressure level at each physical point
    543 
    544      DO l=1,nlayer+1
    545        DO ig=1,ip1jmp1
    546         zplev_new(ig,l) = ap(l)  + bp(l)*ps(ig)
    547        ENDDO
    548      ENDDO
    549 
    550 ! Compute the new values of tracer given the change of pressure and pressure level
    551 
    552      DO nnq=1,nqtot      ! For each tracer
    553       if (nnq.NE.1) then ! If it is not co2
    554         DO l=1,llm-1     ! For each level
    555           DO i=1,ip1jmp1 ! Each point
    556              q(i,l,nnq)=q(i,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))
    557           ENDDO          ! i loop
    558           q(:,llm,nnq)=q(:,llm-1,nnq)
    559         ENDDO            ! l loop
    560       else               ! If it is co2
    561         DO l=1,llm-1     ! For each level
    562           DO i=1,ip1jmp1 ! Each point
    563               q(i,l,nnq)=q(i,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))  &
    564                               + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
    565                                     (zplev_old(ig,l)-zplev_old(ig,l+1))     )  / &
    566                                     (zplev_new(ig,l)-zplev_new(ig,l+1))
    567           ENDDO          ! i loop
    568          q(:,llm,nnq)=q(:,llm-1,nnq)
    569         ENDDO            ! l loop
    570       endif
    571      ENDDO               ! nnq loop
    572 
    573 
    574 !We apply the tendency on the ice
    575 
    576      call evol_h2o_ice_s_slope(qsurf_slope(:,6,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
    577 
    578      call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
    579 
    580 ! FLOW : Make the co2 ice flow
     1008!    II_a update pressure, ice and tracers
     1009!    II_b CO2 glaciers flows
     1010
     1011!------------------------
     1012
     1013   
     1014
    5811015
    5821016       DO ig = 1,ngrid
    … …  
    6191053             flag_co2flow_mesh(ig) = 1.
    6201054            ENDIF ! co2ice > hmax
    621           ENDIF ! iflat
     1055          ENDIF ! iflattsoil_lope
    6221056        ENDDO !islope
    6231057       ENDDO !ig
    6241058
    625 
    626 ! End of FLOW
    627 
    628 ! We recompute the tendencies of co2 given the new surface pressure
    629 
    630     call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,vmr_co2_gcm_phys,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
    631 
     1059      print *, "6 is okay"
     1060
     1061!------------------------
     1062
     1063! II  Run
     1064!    II_a update pressure, ice and tracers
     1065!    II_b CO2 glaciers flows
     1066!    II_c Update surface and soil temperatures
     1067
     1068!------------------------
     1069
     1070
     1071! II_c.1 Update Tsurf
     1072      bool_sublim=0
     1073      Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:)
     1074       DO ig = 1,ngrid
     1075        DO islope = 1,nslope
     1076           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
     1077
     1078
     1079              if(latitude_deg(ig).gt.0) then
     1080                do ig_loop=ig,ngrid
     1081                  DO islope_loop=islope,iflat,-1
     1082                     if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
     1083                            tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
     1084                            bool_sublim=1
     1085                            exit
     1086                     endif
     1087                  enddo
     1088                  if (bool_sublim.eq.1) then
     1089                   exit
     1090                  endif
     1091                enddo
     1092              else
     1093               do ig_loop=ig,1,-1
     1094                  DO islope_loop=islope,iflat
     1095                     if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
     1096                            tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
     1097                            bool_sublim=1
     1098                            exit
     1099                     endif
     1100                  enddo
     1101                  if (bool_sublim.eq.1) then
     1102                   exit
     1103                  endif
     1104                enddo
     1105
     1106              endif
     1107              initial_co2_ice_slope(ig,islope)=0
     1108
     1109              if ((co2ice_slope(ig,islope).lt.1e-5).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then
     1110                 
     1111               albedo_slope(ig,1,islope) = albedo_h2o_frost
     1112               albedo_slope(ig,2,islope) = albedo_h2o_frost
     1113
     1114              else
     1115                           
     1116               albedo_slope(ig,1,islope) = albedodat(ig)
     1117               albedo_slope(ig,2,islope) = albedodat(ig)
     1118               emiss_slope(ig,islope) = emissiv         
     1119
     1120              endif
     1121
     1122           elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2
     1123 
     1124             ave=0.
     1125             do t=1,timelen
     1126              if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then
     1127                 ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
     1128              else
     1129                 ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
     1130              endif
     1131             enddo
     1132
     1133             tsurf_ave_phys(ig,islope)=ave/timelen
     1134           endif
     1135        enddo
     1136       enddo
     1137
     1138      print *, "7 is okay"
     1139   
     1140     do t = 1,timelen
     1141     tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:))
     1142     enddo
     1143! for the start
     1144     do ig = 1,ngrid
     1145      do islope = 1,nslope
     1146     tsurf_slope(ig,islope) =  tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
     1147     enddo
     1148     enddo
     1149
     1150! II_c.2 Update soil temperature
     1151
     1152allocate(TI_locslope(ngrid,nsoilmx_PEM))
     1153allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
     1154allocate(Tsurf_locslope(ngrid))
     1155allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
     1156allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
     1157
     1158
     1159
     1160
     1161print *,"8 is ok"
     1162
     1163! Soil averaged
     1164  do islope = 1,nslope
     1165     TI_locslope(:,:) = TI_PEM(:,:,islope)
     1166     Tsoil_locslope(:,:) = tsoil_PEM(:,:,islope)
     1167     Tsurf_locslope(:) = tsurf_ave_phys(:,islope)
     1168     alph_locslope(:,:) = alph_PEM(:,:,islope)
     1169     beta_locslope(:,:) = beta_PEM(:,:,islope)
     1170     call soil_pem(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
     1171     call soil_pem(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
     1172     do t = 1,timelen
     1173       tsoil_phys_PEM_timeseries(:,:,islope,t) =  tsoil_phys_PEM_timeseries(:,:,islope,t) + (Tsoil_locslope(:,:)- tsoil_PEM(:,:,islope))
     1174     enddo
     1175     TI_PEM(:,:,islope) = TI_locslope(:,:)
     1176     tsoil_PEM(:,:,islope) = Tsoil_locslope(:,:)
     1177     tsurf_ave_phys(:,islope)  =  Tsurf_locslope(:)
     1178     alph_PEM(:,:,islope) = alph_locslope(:,:)
     1179     beta_PEM(:,:,islope) = beta_locslope(:,:)
     1180  enddo
     1181
     1182
     1183      print *, "9 is okay"
     1184
     1185
     1186!!!!!! Check for the time series
     1187
     1188   
     1189      do ig = 1,ngrid
     1190       do islope = 1,nslope
     1191         do iloop = 1,nsoilmx_PEM
     1192 
     1193         if(isnan(tsoil_PEM(ig,iloop,islope))) then
     1194          write(*,*) "is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope)
     1195           stop
     1196         endif
     1197         enddo
     1198        enddo
     1199      enddo
     1200
     1201
     1202write(*,*) "before ice table, no stop"
     1203
     1204! II_c.3 Update the ice table
     1205     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, &
     1206                           ps_phys_timeseries, ice_depth)
     1207       
     1208      print *, "10 is okay"
     1209! II_c.3 Update the soil thermal properties
     1210      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, &
     1211        cell_area,ice_depth,TI_PEM)
     1212
     1213! II_c.4 Update the mass of the regolith adsorbded
     1214   
     1215 
     1216   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, &
     1217                               co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), q_co2_PEM_phys,q_h2o_PEM_phys,co2_adsorbded_phys,delta_co2_adsorbded)
     1218
     1219
     1220!------------------------
     1221
     1222! II  Run
     1223!    II_a update pressure, ice and tracers
     1224!    II_b CO2 glaciers flows
     1225!    II_c Update surface and soil temperatures
     1226!    II_d Update the tendencies and stopping criterion
     1227
     1228!------------------------
     1229
     1230
     1231      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,&     
     1232                               global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
     1233
     1234          print *, "12 is okay"
    6321235      year_iter=year_iter+dt_pem
    6331236
    634 ! We check if we should stop : the criterion is based on the area of ice that is sublimating
    635 
    636       call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,6,:),STOPPING_water,ngrid,initial_h2o_ice)
    637 
    638       call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_slope,global_ave_press_GCM,global_ave_press_new,nslope)
     1237!We check with we should stop
     1238
     1239
     1240      call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
     1241
     1242
     1243          print *, "13 is okay"
     1244      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)
     1245
     1246 print *, "criterion ok"
     1247
    6391248
    6401249      print *, "STOPPING_water", STOPPING_water, "STOPPING1_water", STOPPING_1_water
    6411250      print *, "STOPPING_co2", STOPPING_co2, "STOPPING1_co2", STOPPING_1_co2
    6421251      print *, "year_iter=", year_iter
    643 
    644       year=year+dt_pem
    645 
    646 ! We check if the orbits parameters have changed too much
    647       call orbit_param(year)
    648 
    6491252      if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2)  then
    6501253        exit
    6511254      endif
    6521255
    653 ! Preparing for next iteration
    6541256      global_ave_press_old=global_ave_press_new
    6551257      zplev_old=zplev_new
    6561258
    6571259
    658       enddo
    659 
    660 ! We put the correct thickness of co2ice given the slope
     1260      enddo  ! big loop of the pem
     1261
     1262
     1263!---------------------------- END RUN PEM ---------------------
     1264
     1265
     1266
     1267
     1268!---------------------------- OUTPUT ---------------------
     1269
     1270
     1271
     1272
     1273!------------------------
     1274
     1275! III Output
     1276!    III_a Recomp tendencies for the start and startfi
     1277
     1278!------------------------
     1279
     1280
     1281! III_a.1 Ice update
    6611282      DO ig = 1,ngrid
    6621283          co2ice(ig) = 0.
    … …  
    6681289      ENDDO ! of DO ig=1,ngrid
    6691290
    670 ! We put the correct thickness of water ice given the slope
    671       DO ig = 1,ngrid
    672           qsurf(ig,6) = 0.
     1291      DO ig = 1,ngrid !!!!! TO DOC GHANGE IF QSURF ?ü!?!?!?
     1292          qsurf(ig,igcm_h2o_ice) = 0.
    6731293             DO islope = 1,nslope
    674                  qsurf(ig,6) = qsurf(ig,6) + qsurf_slope(ig,6,islope) &
     1294                 qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
    6751295                            * subslope_dist(ig,islope) /          &
    6761296                           cos(pi*def_slope_mean(islope)/180.)
    … …  
    6801300
    6811301
     1302
     1303! III_a.2 Tsurf and Tsoil update
     1304
     1305
     1306   call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
     1307   tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
     1308
     1309
     1310      DO ig = 1,ngrid
     1311          tsurf(ig) = 0.
     1312         
     1313             DO islope = 1,nslope
     1314                 tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) &
     1315                            * subslope_dist(ig,islope)     
     1316             ENDDO
     1317      ENDDO ! of DO ig=1,ngrid
     1318
     1319      DO ig = 1,ngrid
     1320         DO iloop = 1,nsoilmx
     1321          tsoil(ig,iloop) = 0.
     1322          inertiesoil(ig,iloop) = 0.
     1323             DO islope = 1,nslope
     1324                  tsoil(ig,iloop) =  tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
     1325                            * subslope_dist(ig,islope)   
     1326                  inertiesoil(ig,iloop) =  inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
     1327                            * subslope_dist(ig,islope)         
     1328             ENDDO
     1329          ENDDO
     1330      ENDDO ! of DO ig=1,ngrid
     1331
     1332! III_a.3 Surface optical properties
     1333
     1334    DO ig = 1,ngrid
     1335       DO l = 1,2
     1336        albedo(ig,l) =0.
     1337        DO islope = 1,nslope
     1338          albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) &
     1339                         *subslope_dist(ig,islope)
     1340        ENDDO
     1341       ENDDO
     1342      ENDDO
     1343
     1344
     1345    DO ig = 1,ngrid
     1346
     1347        emis(ig) =0.
     1348        DO islope = 1,nslope
     1349          emis(ig)= emis(ig)+emiss_slope(ig,islope) &
     1350                         *subslope_dist(ig,islope)
     1351        ENDDO
     1352       ENDDO
     1353     
     1354
     1355
    6821356!------------------------
    6831357
    6841358! III Output
    685 !    III_a Write startfi.nc
    686 
    687 !------------------------
    688 
    689 !----------------------------WRITE restart.nc ---------------------
     1359!    III_a Recomp tendencies for the start and startfi
     1360!    III_b Write start and starfi.nc
     1361
     1362!------------------------
     1363
     1364! III_b.1 WRITE restart.nc
    6901365
    6911366      ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys
    … …  
    6931368      ztime_fin=0.
    6941369
     1370         print *, "RVip1jmp1", ip1jmp1
    6951371     allocate(p(ip1jmp1,nlayer+1))
    696      CALL pression (ip1jmp1,ap,bp,ps,p)
    697      CALL massdair(p,masse)
    698 
    699      day_ini=year
    700 
    701      CALL dynredem0("restart_evol.nc", day_ini, phis)
    702 
    703      CALL dynredem1("restart_evol.nc",   &
     1372         print *, "ngrid", ngrid
     1373         print *, "nlayer", nlayer
     1374          CALL pression (ip1jmp1,ap,bp,ps,p)
     1375     print *, "masse 1", masse(1,:)
     1376          CALL massdair(p,masse)
     1377     print *, "masse 2", masse(1,:)
     1378
     1379     do nnq=1,nqtot
     1380        call gr_fi_dyn(llm,ngridmx,iip1,jjp1,q_phys(:,:,nnq),q(:,:,nnq))
     1381     enddo
     1382
     1383
     1384      CALL dynredem0("restart_evol.nc", day_ini, phis)
     1385
     1386      CALL dynredem1("restart_evol.nc",   &
    7041387                time_0,vcov,ucov,teta,q,masse,ps)
    7051388
    7061389
    707 !----------------------------WRITE restartfi.nc ---------------------
    708 
    709      pday=year
    710 
    711      call physdem0("restartfi_evol.nc",longitude,latitude, &
    712                      nsoilmx,ngrid,nlayer,nq,              &
    713                      ptimestep,pday,0.,cell_area,          &
    714                      albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
    715                      hmons,summit,base,nslope,def_slope,   &
    716                      subslope_dist)
    717 
    718     call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,   &
     1390! III_b.2 WRITE restartfi.nc
     1391
     1392      
     1393             call physdem0("restartfi_evol.nc",longitude,latitude, &
     1394                        nsoilmx,ngrid,nlayer,nq,              &
     1395                        ptimestep,pday,0.,cell_area,          &
     1396                        albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
     1397                        hmons,summit,base,nslope,def_slope,  &
     1398                        subslope_dist)
     1399
     1400
     1401          call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
    7191402                     ptimestep,ztime_fin,                        &
    7201403                     tsurf,tsoil,co2ice,albedo,emis,             &
    7211404                     q2,qsurf,tauscaling,totcloudfrac,wstar,     &
    722                      watercap,nslope,co2ice_slope,               &
     1405                     watercap,inertiesoil,nslope,co2ice_slope,   &
    7231406                     tsurf_slope,tsoil_slope, albedo_slope,      &
    724                      emiss_slope,qsurf_slope,watercap_slope)
    725 
    726 
    727       print *, "RV : So far so good"
    728 
    729 
    730 
    731 
    732       deallocate(longitude)
    733       deallocate(latitude)
    734       deallocate(cell_area)
     1407                     emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys)
     1408
     1409
     1410!------------------------
     1411
     1412! III Output
     1413!    III_a Recomp tendencies for the start and startfi
     1414!    III_b Write start and starfi.nc
     1415!    III_c Write start_pem
     1416
     1417!------------------------
     1418
     1419 
     1420
     1421         call pemdem1("restartfi_PEM.nc",ztime_fin,nsoilmx_PEM,ngrid,nslope , &
     1422                      tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys)
     1423
     1424
     1425 
     1426
     1427      print *, "LL & RV : So far so good"
     1428
     1429
     1430
     1431
    7351432
    7361433

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

    r2779 r2794  
    1 !
    2 ! $Id $
    3 !
    4 SUBROUTINE pemetat0(fichnom,min_h2o_ice_s,min_co2_ice_s,iim_input,jjm_input,nlayer,vmr_co2_gcm,ps_GCM,timelen,min_co2_ice_slope,min_h2o_ice_slope,nslope)
    5 
    6       use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
    7                         nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
    8                         nf90_inquire_dimension,nf90_close
    9 
    10       IMPLICIT NONE
    11 
    12 !=======================================================================
    13 !
    14 ! Read initial confitions file
    15 !
    16 !=======================================================================
    17 
    18   include "dimensions.h"
    19 
    20 !===============================================================================
    21 ! Arguments:
    22 
    23 ! INPUT:
    24   CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
    25 
    26   INTEGER :: iim_input,jjm_input,nlayer,nslope
    27 
    28 ! LOCAL:
    29 
    30   REAL, ALLOCATABLE ::  h2o_ice_s(:,:,:)                       ! h2o_ice_s of the concatenated file
    31   REAL, ALLOCATABLE ::  co2_ice_s(:,:,:)                       ! co2_ice_s of the concatenated file
    32 
    33   REAL, ALLOCATABLE ::  h2o_ice_s_slope(:,:,:,:)               ! h2o_ice_s slope of the concatenated file
    34   REAL, ALLOCATABLE ::  co2_ice_s_slope(:,:,:,:)               ! co2_ice_s slope of the concatenated file
    35 
    36   REAL, ALLOCATABLE ::  q_co2_GCM(:,:,:)                       ! vmr of co2 of the first layer in the GCM year
    37 
    38 ! OUTPUT:
    39 
    40   REAL, INTENT(OUT) ::  min_h2o_ice_s(iim_input+1,jjm_input+1) ! Minimum of h2o_ice_s of the year
    41   REAL, INTENT(OUT) ::  min_co2_ice_s(iim_input+1,jjm_input+1) ! Minimum of co2_ice_s of the year
    42 
    43   REAL, INTENT(OUT) ::  min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
    44   REAL, INTENT(OUT) ::  min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
    45 
    46   REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,2676)         ! vmr of co2 of the first layer in the GCM year
    47 
    48   REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,2676)  ! Surface pressure in the GCM year
    49 
    50 !===============================================================================
    51 !   Local Variables
    52   CHARACTER(LEN=256) :: msg, var, modname
    53   INTEGER,PARAMETER :: length=100
    54   INTEGER :: iq, fID, vID, idecal
    55   INTEGER :: ierr
    56   CHARACTER(len=12) :: start_file_type="earth" ! default start file type
    57 
    58   REAL,ALLOCATABLE :: time(:) ! times stored in start
    59   INTEGER :: timelen ! number of times stored in the file
    60   INTEGER :: indextime ! index of selected time
    61 
    62   INTEGER :: edges(4),corner(4)
    63   INTEGER :: i,j,t
    64   real,save :: m_co2, m_noco2, A , B, mmean
    65 
    66   INTEGER :: islope
    67   CHARACTER*2 :: num
    68 
    69 
    70 !-----------------------------------------------------------------------
    71   modname="pemetat0"
    72 
    73       m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
    74       m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
    75       A =(1/m_co2 - 1/m_noco2)
    76       B=1/m_noco2
    77 
    78 !  Open initial state NetCDF file
    79   var=fichnom
    80   CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
    81 
    82       ierr = nf90_inq_varid (fID, "temps", vID)
    83       IF (ierr .NE. nf90_noerr) THEN
    84         write(*,*)"pemetat0: Le champ <temps> est absent"
    85         write(*,*)"pemetat0: J essaie <Time>"
    86         ierr = nf90_inq_varid (fID, "Time", vID)
    87         IF (ierr .NE. nf90_noerr) THEN
    88            write(*,*)"pemetat0: Le champ <Time> est absent"
    89            write(*,*)trim(nf90_strerror(ierr))
    90            CALL ABORT_gcm("pemetat0", "", 1)
    91         ENDIF
    92         ! Get the length of the "Time" dimension
    93         ierr = nf90_inq_dimid(fID,"Time",vID)
    94         ierr = nf90_inquire_dimension(fID,vID,len=timelen)
    95       ELSE   
    96         ! Get the length of the "temps" dimension
    97         ierr = nf90_inq_dimid(fID,"temps",vID)
    98         ierr = nf90_inquire_dimension(fID,vID,len=timelen)
    99       ENDIF
    100 
    101       allocate(co2_ice_s(iim+1,jjm+1,timelen))
    102       allocate(q_co2_GCM(iim+1,jjm+1,timelen))
    103 
    104       allocate(co2_ice_s_slope(iim+1,jjm+1,nslope,timelen))
    105       allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
    106 
    107 
    108 
    109 ! Get all the needed variable of the concatenated file
    110   CALL get_var3("h2o_ice_s"   ,h2o_ice_s)
    111   CALL get_var3("co2ice"   ,co2_ice_s)
    112   CALL get_var3("co2_cropped"   ,q_co2_GCM)
    113   CALL get_var3("ps"   ,ps_GCM)
    114 
    115 ! Slope
     1   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)
     3   
     4   use iostart_PEM, only:  open_startphy, close_startphy, get_field
     5   use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
     6   use comsoil_h, only:  volcapa,inertiedat
     7   use ini_soil_mod, only:  ini_icetable
     8   use soil_evolution_mod, only: soil_pem_CN
     9   use adsorption_mod, only : regolith_co2adsorption
     10  implicit none
     11
     12 
     13  character(len=*), intent(in) :: filename
     14  LOGICAL,intent(in) :: startpem_file
     15  character*8 :: fichnom
     16  integer,intent(in) :: ngrid
     17  integer,intent(in) :: nsoil_GCM
     18  integer,intent(in) :: nsoil_PEM
     19  integer,intent(in) :: nslope
     20  integer,intent(in) :: timelen
     21  real, intent(in) :: tsurf_ave_yr1(ngrid,nslope)                          ! surface temperature at the first year of GCM call
     22  real,intent(in) :: tsurf_ave(ngrid,nslope)                  ! surface temperature at the current year
     23  real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
     24  real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
     25  real,intent(in) :: ps_yr1_inst(ngrid,timelen)
     26  real,intent(in) :: ps_inst(ngrid,timelen)                        ! surface pressure [Pa]
     27  real,intent(in) :: tsurf_inst(ngrid,nslope,timelen) ! soil (mid-layer) temperature
     28  real,intent(in) :: timestep       ! time step
     29  real,intent(in) :: tend_h2oglaciers(ngrid,nslope)
     30  real,intent(in) :: tend_co2glaciers(ngrid,nslope)
     31  real,intent(in) :: co2ice(ngrid,nslope)
     32  real,intent(in) :: waterice(ngrid,nslope)
     33  real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope)
     34
     35! outputs
     36
     37
     38
     39
     40
     41  real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) thermal inertia (SI)
     42  real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) temperature (K)
     43  real,intent(inout) :: ice_table(ngrid,nslope) ! (m)
     44  real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen) ! instantaneous soil (mid-layer) temperature (k)
     45  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed (kg/m^2)
     46  real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption
     47! local
     48   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)
     49   real :: TI_startPEM(ngrid,nsoil_PEM,nslope)
     50   LOGICAL :: found
     51   integer :: iloop,ig,islope,it,isoil
     52   REAL :: TI_breccia = 750.
     53   REAL :: TI_bedrock = 2300.
     54   real :: kcond
     55   real :: delta
     56   CHARACTER*2 :: num
     57   real :: tsoil_saved(ngrid,nsoil_PEM)
     58   real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope)
     59   real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)
     60   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)
     66!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     67!!!
     68!!! Purpose: read start_pem. Need a specific iostart_PEM
     69!!!
     70!!! Order: 1. Thermal Inertia
     71!!!        2. Soil Temperature
     72!!!        3. Ice table
     73!!!        4. Mass of CO2 adsorbed
     74!!!
     75!!! /!\ This order must be respected !
     76!!! Author: LL
     77!!!
     78!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     79
     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
     85!1. Run
     86
     87
     88
     89
     90if (startpem_file) then
     91   ! open pem initial state file:
     92   call open_startphy(filename)
     93 
     94!1. Thermal Inertia
     95! a. General case
    11696
    11797DO islope=1,nslope
    11898  write(num,fmt='(i2.2)') islope
    119   call get_var3("co2ice_slope"//num,co2_ice_s_slope(:,:,islope,:))
    120 ENDDO
     99  call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
     100   if(.not.found) then
     101      write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
     102      write(*,*)'will reconstruct the values of TI_PEM'
     103
     104      do ig = 1,ngrid
     105
     106          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then
     107!!! transition
     108             delta = 50.
     109             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
     110            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
     111                      ((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
     124!! 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))))
     129          do iloop=n_1km+2,nsoil_PEM
     130            TI_PEM(ig,iloop,islope) = TI_bedrock
     131         enddo   
     132      enddo
     133
     134   else
     135     do iloop = nsoil_GCM+1,nsoil_PEM
     136       TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
     137     enddo
     138   endif ! not found
     139ENDDO ! islope
     140
     141print *,'PEMETAT0: THERMAL INERTIA DONE'
     142
     143
     144
     145! b. Special case for inertiedat, inertiedat_PEM
     146call get_field("inertiedat_PEM",inertiedat_PEM,found)
     147if(.not.found) then
     148 do iloop = 1,nsoil_GCM
     149   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
     150 enddo
     151
     152 
     153!!! zone de transition
     154delta = 50.
     155do ig = 1,ngrid
     156if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then
     157inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
     158            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
     159                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
     160
     161
     162 do iloop = nsoil_GCM+2,n_1km
     163       inertiedat_PEM(ig,iloop) = TI_breccia
     164  enddo
     165
     166else
     167   do iloop=nsoil_GCM+1,n_1km
     168      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
     169   enddo
     170endif ! comparison ti breccia
     171enddo!ig
     172
     173!!! zone de transition
     174delta = 1000.
     175do ig = 1,ngrid
     176inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
     177            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
     178                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
     179enddo
     180
     181
     182
     183 do iloop = n_1km+2, nsoil_PEM
     184   do ig = 1,ngrid
     185      inertiedat_PEM(ig,iloop) = TI_bedrock
     186   enddo
     187 enddo
     188endif ! not found
     189
     190
     191!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     192
     193!2. Soil Temperature
    121194
    122195DO islope=1,nslope
    123196  write(num,fmt='(i2.2)') islope
    124   call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_slope(:,:,islope,:))
     197   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
     198  if(.not.found) then
     199      write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
     200      write(*,*)'will reconstruct the values of Tsoil'
     201      do ig = 1,ngrid
     202        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
     203        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))   
     204
     205       do iloop=nsoil_GCM+2,n_1km
     206            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
     207            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
     208      enddo
     209        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
     210        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1)) 
     211
     212       do iloop=n_1km+2,nsoil_PEM
     213            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
     214            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
     215      enddo
     216    enddo
     217
     218   else
     219! predictor corrector: restart from year 1 of the GCM and build the evolution of
     220! tsoil at depth
     221
     222    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)
     225    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)
     227
     228
     229     do iloop = nsoil_GCM+1,nsoil_PEM
     230       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
     231     enddo
     232   endif
     233
     234    do it = 1,timelen
     235        do isoil = nsoil_GCM+1,nsoil_PEM
     236        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
     237        enddo
     238     enddo
     239
     240
     241
     242   
    125243ENDDO
    126 
    127 ! Compute the minimum over the year for each point
    128   min_h2o_ice_s(:,:)=minval(h2o_ice_s,3)
    129   min_co2_ice_s(:,:)=minval(co2_ice_s,3)
    130 
    131   min_co2_ice_slope(:,:,:)=minval(co2_ice_s_slope,4)
    132   min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
    133 
    134 ! By definition, a density is positive, we get rid of the negative values
    135   DO i=1,iim+1
    136     DO j = 1, jjm+1
    137        if (min_co2_ice_s(i,j).LT.0) then
    138           min_h2o_ice_s(i,j)  = 0.
    139           min_co2_ice_s(i,j)  = 0.
    140        endif
    141        DO islope=1,nslope
    142           if (min_co2_ice_slope(i,j,islope).LT.0) then
    143             min_co2_ice_slope(i,j,islope)  = 0.
    144           endif
    145           if (min_h2o_ice_slope(i,j,islope).LT.0) then
    146             min_h2o_ice_slope(i,j,islope)  = 0.
    147           endif
    148        ENDDO
    149     ENDDO
    150   ENDDO
    151 
    152 ! Compute the volume mixing ratio of co2 in the first layer
    153 
    154   DO i=1,iim+1
    155     DO j = 1, jjm+1
    156       DO t = 1, timelen
    157          if (q_co2_GCM(i,j,t).LT.0) then
    158               q_co2_GCM(i,j,t)=1E-10
    159          elseif (q_co2_GCM(i,j,t).GT.1) then
    160               q_co2_GCM(i,j,t)=1.
     244print *,'PEMETAT0: SOIL TEMP  DONE'
     245
     246!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     247
     248!3. Ice Table
     249
     250 
     251   call get_field("ice_table",ice_table,found)
     252   if(.not.found) then
     253      write(*,*)'PEM settings: failed loading <Ice Table>'
     254      write(*,*)'will reconstruct the values of ice table'
     255
     256
     257      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
     259
     260   else
     261   ! update ice table
     262     call computeice_table(timelen,ngrid,nslope,nsoil_PEM,tsoil_inst,tsurf_inst,q_co2,q_h2o, ps_inst, ice_table)
     263
     264
     265   endif
     266
     267print *,'PEMETAT0: ICE TABLE  DONE'
     268
     269
     270
     271!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     272
     273
     274!4. CO2 Adsorption
     275DO islope=1,nslope
     276  write(num,fmt='(i2.2)') islope
     277   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
     278   if(.not.found) then
     279      write(*,*)'PEM settings: failed loading <m_co2_regolith_phys>'
     280      write(*,*)'will reconstruct the values of co2 adsorbded'
     281   m_co2_regolith_phys(:,:,:) = 0.
     282   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
     283
     284   deltam_co2_regolith_phys(:) = 0.
     285   exit
     286   endif
     287ENDDO
     288
     289if (found) then
     290   DO iloop = 1,nsoil_GCM
     291      tsoil_tmp_yr1(:,iloop,:) = tsoil_PEM_yr1(:,iloop,:)
     292
     293   ENDDO 
     294
     295   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
     296
     297
     298endif
     299print *,'PEMETAT0: CO2 adsorption done  '
     300
     301
     302
     303call close_startphy
     304
     305
     306else !No startfi, let's build all by hand
     307
     308
     309
     310!a) Thermal inertia
     311   do islope = 1,nslope
     312      do ig = 1,ngrid
     313
     314          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then
     315!!! transition
     316             delta = 50.
     317             TI_PEM(ig,nsoil_GCM+1,islope) =sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
     318            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
     319                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
     320
     321          do iloop=nsoil_GCM+2,n_1km
     322            TI_PEM(ig,iloop,islope) = TI_breccia
     323         enddo     
     324         else ! we keep the high ti values
     325           do iloop=nsoil_GCM+1,n_1km
     326                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
     327           enddo
     328
    161329         endif
    162          mmean=1/(A*q_co2_GCM(i,j,t) +B)
    163          vmr_co2_gcm(i,j,t) = q_co2_GCM(i,j,t)*mmean/m_co2
     330
     331
     332!! transition
     333             delta = 1000.
     334             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
     335            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
     336                      ((layer_PEM(n_1km+1)-delta)/(TI_breccia**2))))
     337          do iloop=n_1km+2,nsoil_PEM
     338            TI_PEM(ig,iloop,islope) = TI_bedrock
     339         enddo   
     340      enddo
     341
     342enddo
     343
     344
     345 do iloop = 1,nsoil_GCM
     346   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
     347 enddo
     348!!! zone de transition
     349delta = 50.
     350do ig = 1,ngrid
     351      if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then
     352inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
     353            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
     354                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
     355
     356
     357 do iloop = nsoil_GCM+2,n_1km
     358
     359       inertiedat_PEM(ig,iloop) = TI_breccia
     360   
     361 enddo
     362else
     363   do iloop = nsoil_GCM+1,n_1km
     364       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
     365    enddo
     366
     367endif
     368enddo
     369
     370
     371!!! zone de transition
     372delta = 1000.
     373do ig = 1,ngrid
     374inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
     375            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
     376                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
     377enddo
     378
     379
     380
     381 do iloop = n_1km+2, nsoil_PEM
     382   do ig = 1,ngrid
     383      inertiedat_PEM(ig,iloop) = TI_bedrock
     384   enddo
     385 enddo
     386
     387print *,'PEMETAT0: TI  DONE'
     388
     389!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     390
     391
     392!b) Soil temperature
     393    do islope = 1,nslope
     394
     395     write(*,*) "islope=",islope
     396     do ig = 1,ngrid
     397        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
     398        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))
     399
     400       do iloop=nsoil_GCM+2,n_1km
     401            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
     402            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
     403      enddo
     404        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
     405        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))
     406
     407       do iloop=n_1km+2,nsoil_PEM
     408            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
     409            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
     410      enddo
     411!     write(*,*) "ig, islope, T=", ig,islope,tsoil_PEM(ig,:,islope)
     412     enddo
     413   
     414    do it = 1,timelen
     415        do isoil = nsoil_GCM+1,nsoil_PEM
     416        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
     417        enddo
     418     enddo
     419enddo
     420print *,'PEMETAT0: TSOIL DONE  '
     421!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     422!c) Ice table   
     423      do islope = 1,nslope
     424      call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:), & 
     425       tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
     426      enddo
     427!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     428
     429print *,'PEMETAT0: Ice table DONE  '
     430
     431
     432
     433!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     434!d) Regolith adsorbed
     435
     436   m_co2_regolith_phys(:,:,:) = 0.
     437   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
     438   deltam_co2_regolith_phys(:) = 0.
     439
     440print *,'PEMETAT0: CO2 adsorption done  '
     441
     442
     443
     444
     445
     446
     447
     448
     449
     450
     451endif ! of if (startphy_file)
     452
     453!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     454    DO ig = 1,ngrid
     455      DO islope = 1,nslope
     456
     457     write(*,*) 'ig,islope ,ice table=',ig,islope,ice_table(ig,islope)
     458
    164459      ENDDO
    165460    ENDDO
    166   ENDDO
    167 
    168 
    169 
    170 
    171 
    172 
    173   CONTAINS
    174 
    175 SUBROUTINE check_dim(n1,n2,str1,str2)
    176   INTEGER,          INTENT(IN) :: n1, n2
    177   CHARACTER(LEN=*), INTENT(IN) :: str1, str2
    178   CHARACTER(LEN=256) :: s1, s2
    179   IF(n1/=n2) THEN
    180     s1='value of '//TRIM(str1)//' ='
    181     s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
    182     WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
    183     CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
    184   END IF
    185 END SUBROUTINE check_dim
    186 
    187 
    188 SUBROUTINE get_var1(var,v)
    189   CHARACTER(LEN=*), INTENT(IN)  :: var
    190   REAL,             INTENT(OUT) :: v(:)
    191   CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
    192   CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
    193 END SUBROUTINE get_var1
    194 
    195 
    196 SUBROUTINE get_var3(var,v) ! on U grid
    197   CHARACTER(LEN=*), INTENT(IN)  :: var
    198   REAL,             INTENT(OUT) :: v(:,:,:)
    199   CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
    200   CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
    201 
    202 END SUBROUTINE get_var3
    203 
    204 SUBROUTINE get_var4(var,v)
    205   CHARACTER(LEN=*), INTENT(IN)  :: var
    206   REAL,             INTENT(OUT) :: v(:,:,:,:)
    207   CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
    208   CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
    209 END SUBROUTINE get_var4
    210 
    211 SUBROUTINE err(ierr,typ,nam)
    212   INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
    213   CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
    214   CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
    215   IF(ierr==NF90_NoERR) RETURN
    216   SELECT CASE(typ)
    217     CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
    218     CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
    219     CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
    220     CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
    221   END SELECT
    222   CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
    223 END SUBROUTINE err
    224 
    225 END SUBROUTINE pemetat0
     461
     462
     463
     464!! small test
     465    DO ig = 1,ngrid
     466      DO islope = 1,nslope
     467        DO iloop = 1,nsoil_PEM
     468           if(isnan(tsoil_PEM(ig,iloop,islope))) then
     469         write(*,*) "failed nan construction", ig, iloop, islope
     470           stop
     471           endif
     472        ENDDO
     473      ENDDO
     474     ENDDO
     475
     476      write(*,*) "construction ok, no nan"
     477     
     478
     479   END SUBROUTINE

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

    r2779 r2794  
    22! $Id $
    33!
    4 SUBROUTINE recomp_tend_co2_slope(tendencies_co2_ice_phys,vmr_co2_gcm,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
     4SUBROUTINE recomp_tend_co2_slope(tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
    55
    66      IMPLICIT NONE
    … …  
    1616
    1717!   INPUT
    18   INTEGER, intent(in) :: timelen,ngrid,nslope                    ! # of timepoint if diagfi, # grid point, # subslope
    19   REAL, INTENT(in) ::  vmr_co2_gcm(ngrid,timelen)                ! physical point x timelen field : Volume mixing ratio of co2 in the first layer
    20   REAL, intent(in) :: ps_GCM_2(ngrid,timelen)                    ! physical point x timelen field : Surface pressure in the GCM
    21   REAL, intent(in) :: global_ave_press_GCM                       ! Average pressure in the GCM output
    22   REAL, intent(in) :: global_ave_press_new                       ! Actual average pressure
     18  INTEGER, intent(in) :: timelen,ngrid,nslope
     19  REAL, INTENT(in) ::  vmr_co2_gcm(ngrid,timelen)                ! physical point field : Volume mixing ratio of co2 in the first layer
     20  REAL, INTENT(in) ::  vmr_co2_pem(ngrid,timelen)                ! physical point field : Volume mixing ratio of co2 in the first layer
     21  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
     24  REAL, intent(in) :: global_ave_press_GCM
     25  REAL, intent(in) :: global_ave_press_new
     26  REAL, intent(in) ::  tendencies_co2_ice_phys_ini(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    2327
    2428
    2529!   OUTPUT
    26   REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope)  ! physical point field : Evolution of perenial ice over one year
     30  REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    2731
    2832
    … …  
    3034!   ----
    3135
    32   INTEGER :: i,t,islope                                          ! Loop variables
    33   REAL :: eps, sigma, L, beta, alpha, coef, ave                  ! Constant
     36  INTEGER :: i,t,islope
     37  REAL :: eps, sigma, L, beta, alpha, coef, ave
    3438
    3539  eps=0.95
    … …  
    4145  coef=669*24*3600*eps*sigma/L
    4246
     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
    4351! Evolution of the water ice for each physical point
    4452  do i=1,ngrid
    4553    do islope=1,nslope
    4654       ave=0.
    47        if (tendencies_co2_ice_phys(i,islope).LT. 1E-5 .and. tendencies_co2_ice_phys(i,islope).GT.-1E-5) then
    48        else
    49          do t=1,timelen
    50            ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100)))**4  &
    51               -(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)*global_ave_press_GCM/global_ave_press_new/100)))**4
    52          enddo
    53          tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys(i,islope)+coef*ave/timelen
    54        endif
     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  &
     59              -(beta/(alpha-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*global_ave_press_GCM/global_ave_press_new/100)))**4
    5560    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
    5672  enddo
    5773
Note: See TracChangeset for help on using the changeset viewer.