Ignore:
Timestamp:
Feb 16, 2023, 5:29:48 PM (23 months ago)
Author:
romain.vande
Message:

Mars PEM:
Deep cleaning of variables name and allocate.
All the "dyn to phys" grid change is done in subroutines and not in the main program.

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

Legend:

Unmodified
Added
Removed

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

    r2835 r2897  
    22! $Id $
    33!
    4 SUBROUTINE compute_tendencies_slope(tendencies_h2o_ice,min_h2o_ice_Y1,&
    5      min_h2o_ice_Y2,iim_input,jjm_input,ngrid,tendencies_h2o_ice_phys,nslope)
     4SUBROUTINE compute_tendencies_slope(ngrid,nslope,min_ice_Y1,&
     5     min_ice_Y2,tendencies_ice)
    66
    77      IMPLICIT NONE
    … …  
    1818!   INPUT
    1919
    20      INTEGER, intent(in) :: iim_input,jjm_input,ngrid  ,nslope                           ! # of grid points along longitude/latitude/ total
    21      REAL, intent(in) , dimension(iim_input+1,jjm_input+1,nslope):: min_h2o_ice_Y1       ! LON x LAT field : minimum of water ice at each point for the first year
    22      REAL, intent(in) , dimension(iim_input+1,jjm_input+1,nslope):: min_h2o_ice_Y2       ! LON x LAT field : minimum of water ice at each point for the second year
     20     INTEGER, intent(in) :: ngrid, nslope                           ! # of grid points along longitude/latitude/ total
     21     REAL, intent(in) , dimension(ngrid,nslope):: min_ice_Y1       ! LON x LAT field : minimum of water ice at each point for the first year
     22     REAL, intent(in) , dimension(ngrid,nslope):: min_ice_Y2       ! LON x LAT field : minimum of water ice at each point for the second year
    2323
    2424!   OUTPUT
    25      REAL, intent(out) , dimension(iim_input+1,jjm_input+1,nslope) :: tendencies_h2o_ice ! LON x LAT field : difference between the minima = evolution of perenial ice
    26      REAL, intent(out) , dimension(ngrid,nslope)   :: tendencies_h2o_ice_phys            ! physical point field : difference between the minima = evolution of perenial ice
     25     REAL, intent(out) , dimension(ngrid,nslope)   :: tendencies_ice            ! physical point field : difference between the minima = evolution of perenial ice
    2726
    2827!   local:
    2928!   ------
    30 
    31      INTEGER :: i,j,ig0,islope                                                           ! loop variable
     29     INTEGER :: ig,islope                                                           ! loop variable
    3230
    3331!=======================================================================
    … …  
    3533!  We compute the difference
    3634
    37   DO j=1,jjm_input+1
    38     DO i = 1, iim_input
    39        DO islope = 1, nslope
    40          tendencies_h2o_ice(i,j,islope)=min_h2o_ice_Y2(i,j,islope)-min_h2o_ice_Y1(i,j,islope)
    41        enddo
    42     ENDDO
     35  DO ig=1,ngrid
     36    DO islope = 1, nslope
     37      tendencies_ice(ig,islope)=min_ice_Y2(ig,islope)-min_ice_Y1(ig,islope)
     38    enddo
    4339  ENDDO
    4440
    4541!  If the difference is too small; there is no evolution
    46   DO j=1,jjm_input+1
    47     DO i = 1, iim_input
    48        DO islope = 1, nslope
    49          if(abs(tendencies_h2o_ice(i,j,islope)).LT.1.0E-10) then
    50             tendencies_h2o_ice(i,j,islope)=0.
    51          endif
    52        enddo
    53     ENDDO
    54   ENDDO
    55 
    56   DO islope = 1,nslope
    57     CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tendencies_h2o_ice(:,:,islope),tendencies_h2o_ice_phys(:,islope))
     42  DO ig=1,ngrid
     43    DO islope = 1, nslope
     44      if(abs(tendencies_ice(ig,islope)).LT.1.0E-10) then
     45        tendencies_ice(ig,islope)=0.
     46      endif
     47    enddo
    5848  ENDDO
    5949

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

    r2895 r2897  
    133133  REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
    134134  REAL ps(ip1jmp1)                       ! pression  au sol
    135   REAL, dimension(:),allocatable :: ps_phys !(ngrid)                       ! pression  au sol
    136   REAL, dimension(:,:),allocatable :: ps_phys_timeseries !(ngrid x timelen) ! pression  au sol instantannées
    137   REAL, dimension(:,:),allocatable :: ps_phys_timeseries_yr1 !(ngrid x timelen) ! pression  au sol instantannées for the first year of the gcm
     135
     136  REAL, dimension(:),allocatable :: ps_start_GCM !(ngrid)                       ! pression  au sol
     137  REAL, dimension(:,:),allocatable :: ps_timeseries !(ngrid x timelen) ! pression  au sol instantannées
    138138
    139139  REAL masse(ip1jmp1,llm)                ! masse d'air
    … …  
    179179      INTEGER :: criterion_stop  ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
    180180
    181       REAL, dimension(:,:,:),allocatable  :: q_co2_GCM ! Lon x Lat x Time : mass mixing ratio of co2 in the first layer [kg/kg]
    182 
    183181      real,save :: m_co2, m_noco2, A , B, mmean        ! Molar mass of co2, no co2 (Ar, ...), intermediate A, B for computations, mean molar mass of the layer [mol/kg]
    184       real ,allocatable :: vmr_co2_gcm_phys(:,:) ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
     182      real ,allocatable :: vmr_co2_gcm(:,:) ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
    185183      real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times  co2 volume mixing ratio used in the PEM
    186       real ,allocatable :: q_h2o_GCM_phys(:,:)   ! Physics x Times h2o mass mixing ratio in the first layer from the GCM  [kg/kg]
    187       real ,allocatable :: q_co2_GCM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer from the GCM  [kg/kg]
    188       real ,allocatable :: q_co2_PEM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM  [kg/kg]
    189       REAL, ALLOCATABLE ::  ps_GCM(:,:,:)        ! Lon x Lat x Times: surface pressure from the GCM [Pa]
    190       REAL, ALLOCATABLE :: ps_GCM_yr1(:,:,:)     ! Lon x Lat x Times: surface pressure from the 1st year of the GCM [Pa]
    191       REAL, ALLOCATABLE ::  vmr_co2_gcm(:,:,:)   ! Lon x Lat x Times: co2 volumemixing ratio retrieve from the gcm [m^3/m^3]
    192       REAL, ALLOCATABLE ::  q_h2o_GCM(:,:,:)     ! Lon x Lat x Times: h2o volume mixing ratio retrieved from the GCM
    193       REAL ,allocatable ::  q_h2o_PEM_phys(:,:)  ! Physics x Times: h2o mass mixing ratio computed in the PEM
     184      real ,allocatable :: q_co2_PEM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg]
     185      REAL ,allocatable :: q_h2o_PEM_phys(:,:)  ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]
    194186      integer :: timelen                         ! # time samples
    195187      REAL :: ave                                ! intermediate varibale to compute average
    … …  
    204196      REAL ,allocatable :: watercap_slope(:,:)                           ! Physics x Nslope: watercap per slope
    205197      REAL ,allocatable :: watercap_slope_saved                          ! Value saved at the previous time step
    206       REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of co2 ice at each point for the first year [kg/m^2]
    207       REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of co2 ice at each point for the second year [kg/m^2]
    208       REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_1        ! LON x LAT field : minimum of water ice at each point for the first year [kg/m^2]
    209       REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_2        ! LON x LAT field : minimum of water ice at each point for the second year [kg/m^2]
    210       REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope        ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM [kg/m^2]
    211       REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope     ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
     198      REAL , dimension(:,:), allocatable :: min_co2_ice_1                ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
     199      REAL , dimension(:,:), allocatable :: min_co2_ice_2                ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
     200      REAL , dimension(:,:), allocatable :: min_h2o_ice_1                ! ngrid field : minimum of water ice at each point for the first year [kg/m^2]
     201      REAL , dimension(:,:), allocatable :: min_h2o_ice_2                ! ngrid field : minimum of water ice at each point for the second year [kg/m^2]
     202      REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_slope          ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
    212203      REAL, dimension(:,:),allocatable  :: initial_co2_ice_sublim_slope    ! physical point field : Logical array indicating sublimating point of co2 ice
    213204      REAL, dimension(:,:),allocatable  :: initial_h2o_ice_slope           ! physical point field : Logical array indicating if there is water ice at initial state
    214205      REAL, dimension(:,:),allocatable  :: initial_co2_ice_slope           ! physical point field : Logical array indicating if there is co2 ice at initial state
    215       REAL , dimension(:,:,:), allocatable :: tendencies_co2_ice_slope     ! LON x LAT x nslope field : Tendency of evolution of perenial co2 ice over a year
    216       REAL , dimension(:,:,:), allocatable :: tendencies_h2o_ice_slope     ! LON x LAT x slope field : Tendency of evolution of perenial water ice over a year
    217       REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope   ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
    218       REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
    219       REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
     206      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice              ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
     207      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_ini          ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
     208      REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice              ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
    220209      REAL , dimension(:,:), allocatable :: flag_co2flow(:,:)   !(ngrid,nslope)          ! To flag where there is a glacier flow
    221210      REAL , dimension(:), allocatable :: flag_co2flow_mesh(:)  !(ngrid)          ! To flag where there is a glacier flow
    … …  
    223212!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL
    224213
    225      REAL, ALLOCATABLE :: tsurf_ave(:,:,:)          ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
    226      REAL, ALLOCATABLE  :: tsurf_ave_phys(:,:)      ! Physic x LAT x SLOPE field : Averaged Surface Temperature [K]
    227      REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:)        ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
    228      REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:)    ! LON x LAT x SLOPE field : Averaged Soil Temperature  during 1st year of the GCM [K]
    229      REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) ! Physics x SLOPE field : Averaged Soil Temperature  during 1st year [K]
    230      REAL, ALLOCATABLE :: TI_GCM(:,:,:,:)                  ! LON x LAT x SLOPE field : Averaged Thermal Inertia  [SI]
    231      REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:)    ! LON X LAT x SLOPE XTULES field : Surface Temperature in timeseries [K]
    232      REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) ! Physic x SLOPE XTULES field : NOn averaged Surf Temperature [K]
     214     REAL, ALLOCATABLE :: tsurf_ave(:,:)       ! Physic x SLOPE field : Averaged Surface Temperature [K]
     215     REAL, ALLOCATABLE :: tsoil_ave(:,:,:)   ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
     216     REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:) ! Physics x SLOPE field : Averaged Soil Temperature  during 1st year [K]
     217     REAL, ALLOCATABLE :: tsoil_ave_PEM_yr1(:,:,:)
     218     REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:)    ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
    233219
    234220     REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    235      REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:)    !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    236      REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:)     ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
    237      REAL, ALLOCATABLE  :: tsurf_ave_phys_yr1(:,:) ! Physic SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
     221     REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:)    !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
     222     REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
    238223     REAL, ALLOCATABLE :: inertiesoil(:,:)    !Physic x Depth  Thermal inertia of the mesh for restart [SI]
    239224
    240      REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:)  ! Physic x Depth x Slope Averaged GCM Thermal Inertia per slope  [SI]
     225     REAL, ALLOCATABLE :: TI_GCM(:,:,:)       ! Physic x Depth x Slope Averaged GCM Thermal Inertia per slope  [SI]
    241226     REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Same but for the start
    242227
    … …  
    245230     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
    246231     REAL,ALLOCATABLE  :: Tsurf_locslope(:)   ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
    247      REAL,ALLOCATABLE  :: watersurf_density_timeseries(:,:,:,:) ! Lon x Lat x Slope x Times: water surface density, time series [kg/m^3]
    248      REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:,:) ! Lon x Lat x Soil x Slope x Times water soil density, time series [kg /m^3]
    249      REAL,ALLOCATABLE  :: watersurf_density_phys_timeseries(:,:,:) ! Physic  x Slope x Times, water surface density, time series [kg/m^3]
    250      REAL,ALLOCATABLE  :: watersurf_density_phys_ave(:,:)          ! Physic  x Slope, water surface density, yearly averaged [kg/m^3]
    251      REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
    252      REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_ave(:,:,:)          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
     232     REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
     233     REAL,ALLOCATABLE  :: watersurf_density_ave(:,:)          ! Physic  x Slope, water surface density, yearly averaged [kg/m^3]
     234     REAL,ALLOCATABLE  :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
     235     REAL,ALLOCATABLE  :: watersoil_density_PEM_ave(:,:,:)          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
    253236     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)                     ! Surface temperature saved from previous time step [K]
    254237     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)                         ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    … …  
    360343     CALL dynetat0(FILE_NAME_start,vcov,ucov, &
    361344                    teta,q,masse,ps,phis, time_0)
     345
     346     allocate(ps_start_GCM(ngrid))
     347     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
    362348
    363349     CALL iniconst !new
    … …  
    478464     PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
    479465
    480 
    481466     allocate(flag_co2flow(ngrid,nslope))
    482467     allocate(flag_co2flow_mesh(ngrid))
    … …  
    484469       flag_co2flow(:,:) = 0.     
    485470       flag_co2flow_mesh(:) = 0.
    486 
    487471
    488472!---------------------------- READ GCM data ---------------------
    … …  
    492476!    I_b READ of start_evol.nc and starfi_evol.nc
    493477!    I_c Subslope parametrisation
    494 !    I_d READ GCM data and convert to the physical grid
     478!    I_d READ GCM data
    495479
    496480!------------------------
    … …  
    500484     call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
    501485
    502 
    503      allocate(vmr_co2_gcm(iim+1,jjm+1,timelen))
    504      allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
    505      allocate(q_co2_GCM(iim+1,jjm+1,timelen))
    506      allocate(ps_GCM(iim+1,jjm+1,timelen))
    507      allocate(ps_GCM_yr1(iim+1,jjm+1,timelen))
    508      allocate(min_co2_ice_slope_1(iim+1,jjm+1,nslope))
    509      allocate(min_h2o_ice_slope_1(iim+1,jjm+1,nslope))
    510      allocate(tsurf_ave(iim+1,jjm+1,nslope))
    511      allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope))
    512      allocate(tsurf_ave_phys(ngrid,nslope))
    513      allocate(tsurf_ave_phys_yr1(ngrid,nslope))
    514      allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen))
    515      allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen))
    516      allocate(co2_ice_GCM_phys_slope(ngrid,nslope,timelen))
    517      allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen))
     486     allocate(vmr_co2_gcm(ngrid,timelen))
     487     allocate(ps_timeseries(ngrid,timelen))
     488     allocate(min_co2_ice_1(ngrid,nslope))
     489     allocate(min_h2o_ice_1(ngrid,nslope))
     490     allocate(min_co2_ice_2(ngrid,nslope))
     491     allocate(min_h2o_ice_2(ngrid,nslope))
     492     allocate(tsurf_ave_yr1(ngrid,nslope))
     493     allocate(tsurf_ave(ngrid,nslope))
     494     allocate(tsoil_ave_yr1(ngrid,nsoilmx,nslope))
     495     allocate(tsoil_ave(ngrid,nsoilmx,nslope))
     496     allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
     497     allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
     498     allocate(TI_GCM(ngrid,nsoilmx,nslope))
     499     allocate(q_co2_PEM_phys(ngrid,timelen))
     500     allocate(q_h2o_PEM_phys(ngrid,timelen))
     501     allocate(co2_ice_GCM_slope(ngrid,nslope,timelen))
     502     allocate(watersurf_density_ave(ngrid,nslope))
     503     allocate(watersoil_density_timeseries(nslope,nsoilmx,nslope,timelen))
     504
     505     allocate(tsoil_ave_PEM_yr1(ngrid,nsoilmx_PEM,nslope))
    518506     allocate(Tsurfave_before_saved(ngrid,nslope))
    519      allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
    520      allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
    521      allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
    522      allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
    523      allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
    524507     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
     508     allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
     509     allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
    525510     allocate(delta_co2_adsorbded(ngrid))
    526511     allocate(delta_h2o_adsorbded(ngrid))
    527      allocate(watersurf_density_timeseries(iim+1,jjm+1,nslope,timelen))
    528      allocate(watersoil_density_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
    529      allocate(watersurf_density_phys_timeseries(ngrid,nslope,timelen))
    530      allocate(watersurf_density_phys_ave(ngrid,nslope))
    531      allocate(watersoil_density_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    532      allocate(watersoil_density_phys_PEM_ave(ngrid,nsoilmx_PEM,nslope))
     512     allocate(vmr_co2_pem_phys(ngrid,timelen))
    533513
    534514     print *, "Downloading data Y1..."
    535515
    536      call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm,vmr_co2_gcm,ps_GCM_yr1,min_co2_ice_slope_1,min_h2o_ice_slope_1,&      
    537                        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, &     
    538                        watersurf_density_timeseries,watersoil_density_timeseries)
     516     call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1,&   
     517                       tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM_slope, &     
     518                       watersurf_density_ave,watersoil_density_timeseries)
    539519
    540520! 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
    541521
    542522     print *, "Downloading data Y1 done"
    543 
    544      allocate(min_co2_ice_slope_2(iim+1,jjm+1,nslope))
    545      allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
    546 
    547523     print *, "Downloading data Y2"
    548524
    549      call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm,vmr_co2_gcm,ps_GCM,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
    550                   nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &     
    551                   watersurf_density_timeseries,watersoil_density_timeseries)
     525     call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
     526                  tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM_slope, &     
     527                  watersurf_density_ave,watersoil_density_timeseries)
    552528
    553529     print *, "Downloading data Y2 done"
    554 
    555 ! The variables in the dynamic grid are transfered to the physical grid
    556 
    557      allocate(vmr_co2_gcm_phys(ngrid,timelen))
    558      allocate(vmr_co2_pem_phys(ngrid,timelen))
    559      allocate(q_h2o_GCM_phys(ngrid,timelen))
    560      allocate(q_h2o_PEM_phys(ngrid,timelen))
    561      allocate(q_co2_GCM_phys(ngrid,timelen))
    562      allocate(q_co2_PEM_phys(ngrid,timelen))
    563      allocate(ps_phys(ngrid))
    564      allocate(ps_phys_timeseries(ngrid,timelen))
    565      allocate(ps_phys_timeseries_yr1(ngrid,timelen))
    566 
    567      CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys)
    568      CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys)
    569      CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys)
    570      call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
    571      call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
    572      call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
    573      CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
    574      CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
    575 
    576      deallocate(vmr_co2_gcm)
    577      deallocate(q_h2o_GCM)
    578      deallocate(q_co2_GCM)
    579      deallocate(ps_GCM)
    580      deallocate(ps_GCM_yr1)
    581      deallocate(tsurf_ave)
    582      deallocate(tsurf_ave_yr1)
    583 
    584      q_co2_PEM_phys(:,:)=  q_co2_GCM_phys(:,:)
    585      q_h2o_PEM_phys(:,:)=  q_h2o_GCM_phys(:,:)
    586530
    587531!------------------------
    … …  
    602546      call end_adsorption_h_PEM
    603547      call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
     548
    604549      allocate(ice_depth(ngrid,nslope))
    605550      ice_depth(:,:) = 0.
    606       allocate(TI_GCM_phys(ngrid,nsoilmx,nslope))
    607 
    608       DO islope = 1,nslope
    609         if(soil_pem) then
    610           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope))
    611         endif !soil_pem
    612         DO t=1,timelen
    613           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
    614           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
    615         enddo
     551
     552   if(soil_pem) then
     553      call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM,TI_PEM)
     554      DO l=1,nsoilmx
     555        tsoil_ave_PEM_yr1(:,l,:)=tsoil_ave_yr1(:,l,:)
     556        tsoil_PEM(:,l,:)=tsoil_ave(:,l,:)
     557        tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:)
     558        watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,l,:,:)
    616559      ENDDO
    617 
    618       deallocate(co2_ice_GCM_slope)
    619       deallocate(TI_GCM)
    620       deallocate(tsurf_GCM_timeseries)
    621 
    622    if(soil_pem) then
    623       call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
    624       DO islope = 1,nslope
    625         DO t=1,timelen
    626           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersurf_density_timeseries(:,:,islope,t),watersurf_density_phys_timeseries(:,islope,t))
    627         ENDDO
    628         DO l=1,nsoilmx
    629            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
    630            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
    631            DO t=1,timelen
    632              CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
    633            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersoil_density_timeseries(:,:,l,islope,t),watersoil_density_phys_PEM_timeseries(:,l,islope,t))
    634            ENDDO
    635         ENDDO
    636         DO l=nsoilmx+1,nsoilmx_PEM
    637            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
    638            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
    639            DO t=1,timelen
    640            watersoil_density_phys_PEM_timeseries(:,l,islope,t) =  watersoil_density_phys_PEM_timeseries(:,nsoilmx,islope,t)
    641            ENDDO
    642         ENDDO
     560      DO l=nsoilmx+1,nsoilmx_PEM
     561        tsoil_ave_PEM_yr1(:,l,:)=tsoil_ave_yr1(:,nsoilmx,:)
     562        tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:)
     563        watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:)
    643564      ENDDO
    644       watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
    645       watersurf_density_phys_ave(:,:) = SUM(watersurf_density_phys_timeseries(:,:,:),3)/timelen
    646       deallocate(watersurf_density_timeseries)
    647       deallocate(watersurf_density_phys_timeseries)
    648       deallocate(watersoil_density_timeseries)
     565      watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
    649566      deallocate(tsoil_ave_yr1)
    650567      deallocate(tsoil_ave)
    … …  
    664581!----- Compute tendencies from the PCM run
    665582
    666      allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
    667      allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
    668      allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope))
    669      allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope))
    670      allocate(tendencies_h2o_ice_phys_slope(ngrid,nslope))
     583     allocate(tendencies_co2_ice(ngrid,nslope))
     584     allocate(tendencies_co2_ice_ini(ngrid,nslope))
     585     allocate(tendencies_h2o_ice(ngrid,nslope))
    671586
    672587!  Compute the tendencies of the evolution of ice over the years
    673588
    674       call compute_tendencies_slope(tendencies_co2_ice_slope,min_co2_ice_slope_1,&
    675              min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope)
    676 
    677       tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:)
    678 
    679       call compute_tendencies_slope(tendencies_h2o_ice_slope,min_h2o_ice_slope_1,&
    680              min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope)
     589      call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,&
     590             min_co2_ice_2,tendencies_co2_ice)
     591
     592      tendencies_co2_ice_ini(:,:)=tendencies_co2_ice(:,:)
     593
     594      call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,&
     595             min_h2o_ice_2,tendencies_h2o_ice)
     596
     597     deallocate(min_co2_ice_1)
     598     deallocate(min_co2_ice_2)
     599     deallocate(min_h2o_ice_1)
     600     deallocate(min_h2o_ice_2)
    681601
    682602!------------------------
    … …  
    705625    Total_surface=Total_surface+cell_area(i)
    706626    do islope=1,nslope
    707       if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
     627      if (tendencies_co2_ice(i,islope).LT.0) then
    708628         initial_co2_ice_sublim_slope(i,islope)=1.
    709629         ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
    … …  
    716636         initial_co2_ice_slope(i,islope)=0.         
    717637      endif
    718       if (tendencies_h2o_ice_phys_slope(i,islope).LT.0) then
     638      if (tendencies_h2o_ice(i,islope).LT.0) then
    719639         initial_h2o_ice_slope(i,islope)=1.
    720640         ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
    … …  
    733653     DO l=1,nlayer+1
    734654       DO ig=1,ngrid
    735          zplev_gcm(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
     655         zplev_gcm(ig,l) = ap(l)  + bp(l)*ps_start_GCM(ig)
    736656       ENDDO
    737657     ENDDO
    … …  
    739659     global_ave_press_old=0.
    740660     do i=1,ngrid
    741        global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface
     661       global_ave_press_old=global_ave_press_old+cell_area(i)*ps_start_GCM(i)/Total_surface
    742662     enddo
    743663
    … …  
    760680!---------------------------- Read the PEMstart ---------------------
    761681
    762       call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
    763       tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsoil_phys_PEM_timeseries,&
    764       tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,&
    765       watersurf_density_phys_ave,watersoil_density_phys_PEM_ave, &
     682      call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_PEM_yr1,tsoil_PEM,ice_depth, &
     683      tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,&
     684      tendencies_h2o_ice,tendencies_co2_ice,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,&
     685      watersurf_density_ave,watersoil_density_PEM_ave, &
    766686      co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
    767687
    768688    do ig = 1,ngrid
    769689      do islope = 1,nslope
    770 !          qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
     690          qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    771691      enddo
    772692    enddo
    … …  
    791711     write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
    792712     endif ! adsorption
    793      deallocate(tsoil_ave_phys_yr1) 
    794      deallocate(tsurf_ave_phys_yr1)
    795      deallocate(ps_phys_timeseries_yr1) 
     713     deallocate(tsoil_ave_PEM_yr1) 
     714     deallocate(tsurf_ave_yr1)
    796715
    797716!------------------------
    … …  
    836755     do i=1,ngrid
    837756       do islope=1,nslope
    838            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
     757           global_ave_press_new=global_ave_press_new-g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
    839758      enddo
    840759     enddo
    … …  
    856775     DO l=1,nlayer+1
    857776       DO ig=1,ngrid
    858          zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
     777         zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_timeseries(ig,:)
    859778       ENDDO
    860779     ENDDO
    … …  
    862781! II.a.3. Surface pressure timeseries
    863782     print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
    864      do i = 1,ngrid
    865        ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old
     783     do ig = 1,ngrid
     784       ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
    866785     enddo
    867786
    … …  
    871790     do l=1,nlayer+1
    872791       do ig=1,ngrid
    873          zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
     792         zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_timeseries(ig,:)
    874793       enddo
    875794     enddo
    … …  
    916835! II.b. Evolution of the ice
    917836      print *, "Evolution of h2o ice"
    918      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)
    919 
     837     call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
     838
     839     DO islope=1, nslope
     840       write(str2(1:2),'(i2.2)') islope
     841       call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf_slope(:,igcm_h2o_ice,islope))
     842       call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
     843     ENDDO
    920844
    921845      print *, "Evolution of co2 ice"
    922      call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
    923 
    924       DO islope=1, nslope
    925         write(str2(1:2),'(i2.2)') islope
    926         call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf_slope(:,igcm_h2o_ice,islope))
    927         call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice_phys_slope(:,islope))
    928       ENDDO
     846     call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
    929847
    930848!------------------------
    … …  
    939857      print *, "Co2 glacier flows"
    940858
    941     call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_phys_timeseries,&
     859    call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,&
    942860                         global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
    943861
    944       DO islope=1, nslope
    945         write(str2(1:2),'(i2.2)') islope
    946         call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,co2ice_slope(:,islope))
    947         call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice_phys_slope(:,islope))
    948       ENDDO
     862     DO islope=1, nslope
     863       write(str2(1:2),'(i2.2)') islope
     864       call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,co2ice_slope(:,islope))
     865       call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
     866     ENDDO
    949867
    950868!------------------------
    … …  
    962880      print *, "Updating the new Tsurf"
    963881      bool_sublim=0
    964       Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:)
     882      Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
    965883      DO ig = 1,ngrid
    966884        DO islope = 1,nslope
    … …  
    970888                DO islope_loop=islope,iflat,-1
    971889                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
    972                     tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
     890                    tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
    973891                    bool_sublim=1
    974892                    exit
    … …  
    983901                DO islope_loop=islope,iflat
    984902                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
    985                     tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
     903                    tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
    986904                    bool_sublim=1
    987905                    exit
    … …  
    1002920              emiss_slope(ig,islope) = emissiv         
    1003921            endif
    1004           elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice_phys_slope(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
     922          elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
    1005923            ave=0.
    1006924            do t=1,timelen
    1007               if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-3) then
    1008                 ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
     925              if(co2_ice_GCM_slope(ig,islope,t).gt.1e-3) then
     926                ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
    1009927              else
    1010                 ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
     928                ave = ave + tsurf_GCM_timeseries(ig,islope,t)
    1011929              endif
    1012930            enddo
    1013             tsurf_ave_phys(ig,islope)=ave/timelen
     931            tsurf_ave(ig,islope)=ave/timelen
    1014932          endif
    1015933        enddo
    … …  
    1017935
    1018936      do t = 1,timelen
    1019         tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:))
     937        tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) +( tsurf_ave(:,:) -Tsurfave_before_saved(:,:))
    1020938      enddo
    1021939      ! for the start
    1022940      do ig = 1,ngrid
    1023941        do islope = 1,nslope
    1024           tsurf_slope(ig,islope) =  tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
     942          tsurf_slope(ig,islope) =  tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope))
    1025943        enddo
    1026944      enddo
    … …  
    1045963     do t = 1,timelen
    1046964       Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
    1047        Tsurf_locslope(:) =  tsurf_phys_GCM_timeseries(:,islope,t)
     965       Tsurf_locslope(:) =  tsurf_GCM_timeseries(:,islope,t)
    1048966       call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
    1049967       call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
    … …  
    1051969       do ig = 1,ngrid
    1052970         do isoil = 1,nsoilmx_PEM
    1053           watersoil_density_phys_PEM_timeseries(ig,isoil,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
     971          watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
    1054972          if(isnan(Tsoil_locslope(ig,isoil))) then
    1055973            call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1)
    … …  
    1060978  enddo
    1061979     tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
    1062      watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
     980     watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
    1063981
    1064982  print *, "Update of soil temperature done"
    … …  
    1070988
    1071989! II_d.3 Update the ice table
    1072      call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
     990     call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,ice_depth)
    1073991       
    1074992      print *, "Update soil propreties"
    1075993! II_d.4 Update the soil thermal properties
    1076       call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
     994      call update_soil(ngrid,nslope,nsoilmx,nsoilmx_PEM,tendencies_h2o_ice,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
    1077995        ice_depth,TI_PEM)
    1078996
    1079997! II_d.5 Update the mass of the regolith adsorbded
    1080998     if(adsorption_pem) then
    1081     call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),co2ice_slope, &
    1082                              tsoil_PEM,TI_PEM,ps_phys_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
     999       call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf_slope(:,igcm_h2o_ice,:),co2ice_slope, &
     1000                             tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
    10831001                             h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
    1084  
    10851002     endif
    10861003  endif !soil_pem
    … …  
    10981015
    10991016      print *, "Adaptation of the new co2 tendencies to the current pressure"
    1100       call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,co2ice_slope,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,&     
     1017      call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,co2ice_slope,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,&     
    11011018                               global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
    11021019
    … …  
    11521069      endif
    11531070
    1154 
    1155 
    11561071      global_ave_press_old=global_ave_press_new
    11571072
    … …  
    11851100
    11861101! H2O ice
    1187 
    11881102     DO ig=1,ngrid
    11891103       if(watercaptag(ig)) then
    … …  
    12461160   if(soil_pem) then
    12471161     fluxgeo_slope(:,:) = fluxgeo
    1248      call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
     1162     call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM)
    12491163     tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
    12501164   else
    1251       TI_GCM_phys(:,:,:)=TI_GCM_start(:,:,:)
     1165      TI_GCM(:,:,:)=TI_GCM_start(:,:,:)
    12521166   endif !soil_pem
    12531167
    … …  
    12601174            tsoil(ig,iloop) =  tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
    12611175                            * subslope_dist(ig,islope)   
    1262             inertiesoil(ig,iloop) =  inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
     1176            inertiesoil(ig,iloop) =  inertiesoil(ig,iloop) + TI_GCM(ig,iloop,islope) &
    12631177                            * subslope_dist(ig,islope)         
    12641178          ENDDO
    … …  
    12861200    ENDDO
    12871201
    1288 
    12891202      DO ig = 1,ngrid
    12901203        tsurf(ig) = 0.
    … …  
    13041217
    13051218     do i = 1,ngrid
    1306        ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_GCM
     1219       ps_start_GCM(i)=ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM
    13071220     enddo
    13081221
    … …  
    13121225     do l=1,nlayer+1
    13131226       do ig=1,ngrid
    1314          zplev_new(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
     1227         zplev_new(ig,l) = ap(l)  + bp(l)*ps_start_GCM(ig)
    13151228       enddo
    13161229     enddo
    … …  
    13961309                     watercap,inertiesoil,nslope,co2ice_slope,   &
    13971310                     tsurf_slope,tsoil_slope, albedo_slope,      &
    1398                      emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys)
     1311                     emiss_slope,qsurf_slope,watercap_slope, TI_GCM)
    13991312#else
    14001313     call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
    … …  
    14241337        call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
    14251338                      tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
     1339
    14261340        call info_run_PEM(year_iter, criterion_stop)
    14271341
    … …  
    14301344      print *, "LL & RV : So far so good"
    14311345
     1346     deallocate(vmr_co2_gcm)
     1347     deallocate(ps_timeseries)
     1348     deallocate(tsurf_GCM_timeseries)
     1349     deallocate(q_co2_PEM_phys)
     1350     deallocate(q_h2o_PEM_phys)
     1351     deallocate(co2_ice_GCM_slope)
     1352     deallocate(watersurf_density_ave)
     1353     deallocate(watersoil_density_timeseries)
     1354     deallocate(Tsurfave_before_saved)
     1355     deallocate(tsoil_phys_PEM_timeseries)
     1356     deallocate(watersoil_density_PEM_timeseries)
     1357     deallocate(watersoil_density_PEM_ave)
     1358     deallocate(delta_co2_adsorbded)
     1359     deallocate(delta_h2o_adsorbded)
     1360     deallocate(vmr_co2_pem_phys)
     1361
    14321362END PROGRAM pem

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

    r2895 r2897  
    22! $Id $
    33!
    4 SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,vmr_co2_gcm,ps_GCM, &
    5              min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope, &
    6              watersurf_density,watersoil_density)
     4SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries, &
     5             min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2,q_h2o,co2_ice_slope, &
     6             watersurf_density_ave,watersoil_density)
    77
    88      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
    … …  
    2828  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
    2929  INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
    30   INTEGER :: iim_input,jjm_input,nslope            ! number of points in the lat x lon dynamical grid, number of subgrid slopes
    31 
     30  INTEGER :: iim_input,jjm_input,ngrid,nslope            ! number of points in the lat x lon dynamical grid, number of subgrid slopes
    3231! Ouputs
    33   REAL, INTENT(OUT) ::  min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2 ice  per slope of the year [kg/m^2]
    34   REAL, INTENT(OUT) ::  min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of h2o ice per slope of the year [kg/m^2]
    35   REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)      ! CO2 volume mixing ratio in the first layer  [mol/m^3]
    36   REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,timelen)        ! H2O mass mixing ratio in the first layer [kg/m^3]
    37   REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,timelen)        ! CO2 mass mixing ratio in the first layer [kg/m^3]
    38   REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,timelen)          ! Surface Pressure [Pa]
    39   REAL, INTENT(OUT) ::  co2_ice_slope(iim_input+1,jjm_input+1,nslope,timelen) ! co2 ice amount per  slope of the year [kg/m^2]
    40 
     32  REAL, INTENT(OUT) ::  min_co2_ice(ngrid,nslope) ! Minimum of co2 ice  per slope of the year [kg/m^2]
     33  REAL, INTENT(OUT) ::  min_h2o_ice(ngrid,nslope) ! Minimum of h2o ice per slope of the year [kg/m^2]
     34  REAL, INTENT(OUT)  :: vmr_co2_gcm_phys(ngrid,timelen) ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
     35  REAL, INTENT(OUT) ::  ps_timeseries(ngrid,timelen)! Surface Pressure [Pa]
     36  REAL, INTENT(OUT) ::  q_co2(ngrid,timelen)        ! CO2 mass mixing ratio in the first layer [kg/m^3]
     37  REAL, INTENT(OUT) ::  q_h2o(ngrid,timelen)        ! H2O mass mixing ratio in the first layer [kg/m^3]
     38  REAL, INTENT(OUT) ::  co2_ice_slope(ngrid,nslope,timelen) ! co2 ice amount per  slope of the year [kg/m^2]
    4139!SOIL
    42   REAL, INTENT(OUT) ::  tsurf_ave(iim_input+1,jjm_input+1,nslope)         ! Average surface temperature of the concatenated file [K]
    43   REAL, INTENT(OUT) ::  tsoil_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average soil temperature of the concatenated file [K]
    44   REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen)                  ! Surface temperature of the concatenated file, time series [K]
    45   REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)         ! Soil temperature of the concatenated file, time series [K]
    46   REAL , INTENT(OUT) ::  watersurf_density(iim_input+1,jjm_input+1,nslope,timelen)         ! Water density at the surface, time series [kg/m^3]
    47   REAL , INTENT(OUT) ::  watersoil_density(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
    48   REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)                          ! Thermal Inertia  of the concatenated file, times series [SI]
    49   REAL, INTENT(OUT) ::  TI_ave(iim_input+1,jjm_input+1,nsoilmx,nslope)                     ! Average Thermal Inertia  of the concatenated file [SI]
     40  REAL, INTENT(OUT) ::  tsurf_ave(ngrid,nslope)         ! Average surface temperature of the concatenated file [K]
     41  REAL, INTENT(OUT) ::  tsoil_ave(ngrid,nsoilmx,nslope) ! Average soil temperature of the concatenated file [K]
     42  REAL ,INTENT(OUT) ::  tsurf_gcm(ngrid,nslope,timelen)                  ! Surface temperature of the concatenated file, time series [K]
     43  REAL , INTENT(OUT) ::  tsoil_gcm(ngrid,nsoilmx,nslope,timelen)         ! Soil temperature of the concatenated file, time series [K]
     44  REAL , INTENT(OUT) ::  watersurf_density_ave(ngrid,nslope)             ! Water density at the surface [kg/m^3]
     45  REAL , INTENT(OUT) ::  watersoil_density(ngrid,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
     46  REAL, INTENT(OUT) ::  TI_ave(ngrid,nsoilmx,nslope)                     ! Average Thermal Inertia  of the concatenated file [SI]
    5047!===============================================================================
    5148!   Local Variables
    … …  
    5956
    6057  INTEGER :: edges(4),corner(4)
    61   INTEGER :: i,j,t                                                     ! loop variables
     58  INTEGER :: i,j,l,t                                                     ! loop variables
    6259  real,save :: m_co2, m_noco2, A , B, mmean                            ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
    6360
    6461  INTEGER :: islope                                                    ! loop for variables
    6562  CHARACTER*2 :: num                                                   ! for reading sloped variables
    66   REAL, ALLOCATABLE ::  h2o_ice_s(:,:,:)                               ! h2o ice, mesh averaged, of the concatenated file [kg/m^2]
    67   REAL, ALLOCATABLE ::  co2_ice_s(:,:,:)                               ! co2 ice, mesh averaged, of the concatenated file [kg/m^2]
    68   REAL, ALLOCATABLE ::  h2o_ice_s_slope(:,:,:,:)                       ! h2o ice per slope of the concatenated file [kg/m^2]
    69   REAL, ALLOCATABLE ::  watercap_slope(:,:,:,:)
     63  REAL ::  h2o_ice_s_dyn(iim_input+1,jjm_input+1,nslope,timelen)       ! h2o ice per slope of the concatenated file [kg/m^2]
     64  REAL ::  watercap_slope(iim_input+1,jjm_input+1,nslope,timelen)
     65  REAL ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)                ! CO2 volume mixing ratio in the first layer  [mol/m^3]
     66  REAL ::  ps_GCM(iim_input+1,jjm_input+1,timelen)                     ! Surface Pressure [Pa]
     67  REAL ::  min_co2_ice_dyn(iim_input+1,jjm_input+1,nslope)
     68  REAL ::  min_h2o_ice_dyn(iim_input+1,jjm_input+1,nslope)
     69  REAL ::  tsurf_ave_dyn(iim_input+1,jjm_input+1,nslope)               ! Average surface temperature of the concatenated file [K]
     70  REAL ::  tsoil_ave_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope)       ! Average soil temperature of the concatenated file [K]
     71  REAL ::  tsurf_gcm_dyn(iim_input+1,jjm_input+1,nslope,timelen)       ! Surface temperature of the concatenated file, time series [K]
     72  REAL ::  tsoil_gcm_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)! Soil temperature of the concatenated file, time series [K]
     73  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)      ! Thermal Inertia  of the concatenated file, times series [SI]
     74  REAL ::  TI_ave_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope)          ! Average Thermal Inertia  of the concatenated file [SI]
     75  REAL ::  q_co2_dyn(iim_input+1,jjm_input+1,timelen)                  ! CO2 mass mixing ratio in the first layer [kg/m^3]
     76  REAL ::  q_h2o_dyn(iim_input+1,jjm_input+1,timelen)                  ! H2O mass mixing ratio in the first layer [kg/m^3]
     77  REAL ::  co2_ice_slope_dyn(iim_input+1,jjm_input+1,nslope,timelen)  ! co2 ice amount per  slope of the year [kg/m^2]
     78  REAL ::  watersurf_density_dyn(iim_input+1,jjm_input+1,nslope,timelen)! Water density at the surface, time series [kg/m^3]
     79  REAL ::  watersurf_density(ngrid,nslope,timelen)                     ! Water density at the surface, time series [kg/m^3]
     80  REAL ::  watersoil_density_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
     81
    7082!-----------------------------------------------------------------------
    7183  modname="read_data_gcm"
    … …  
    7688      B=1/m_noco2
    7789
    78       allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
    79 
    8090  print *, "Opening ", fichnom, "..."
    8191
    … …  
    8696     print *, "Downloading data for vmr co2..."
    8797
    88   CALL get_var3("co2_cropped"   ,q_co2_GCM)
     98  CALL get_var3("co2_cropped"   ,q_co2_dyn)
    8999
    90100     print *, "Downloading data for vmr co2 done"
    91101     print *, "Downloading data for vmr h20..."
    92102
    93   CALL get_var3("h2o_cropped"   ,q_h2o_GCM)
     103  CALL get_var3("h2o_cropped"   ,q_h2o_dyn)
    94104
    95105     print *, "Downloading data for vmr h2o done"
    … …  
    106116DO islope=1,nslope
    107117  write(num,fmt='(i2.2)') islope
    108   call get_var3("co2ice_slope"//num,co2_ice_slope(:,:,islope,:))
     118  call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
    109119ENDDO
    110120
    … …  
    114124DO islope=1,nslope
    115125  write(num,fmt='(i2.2)') islope
    116   call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_slope(:,:,islope,:))
     126  call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
    117127ENDDO
    118128
    … …  
    122132DO islope=1,nslope
    123133       write(num,fmt='(i2.2)') islope
    124 !       call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
    125         watercap_slope(:,:,:,:)= 0.
     134       call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
     135!        watercap_slope(:,:,:,:)= 0.
    126136ENDDO           
    127137     print *, "Downloading data for watercap_slope done"
    … …  
    131141DO islope=1,nslope
    132142  write(num,fmt='(i2.2)') islope
    133   call get_var3("tsurf_slope"//num,tsurf_gcm(:,:,islope,:))
     143  call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
    134144ENDDO
    135145
    … …  
    142152DO islope=1,nslope
    143153  write(num,fmt='(i2.2)') islope
    144   call get_var4("tsoil_slope"//num,tsoil_gcm(:,:,:,islope,:))
     154  call get_var4("tsoil_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
    145155ENDDO
    146156
    … …  
    159169DO islope=1,nslope
    160170  write(num,fmt='(i2.2)') islope
    161   call get_var4("Waterdensity_soil_slope"//num,watersoil_density(:,:,:,islope,:))
     171  call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
    162172ENDDO
    163173
    … …  
    168178DO islope=1,nslope
    169179  write(num,fmt='(i2.2)') islope
    170   call get_var3("Waterdensity_surface"//num,watersurf_density(:,:,islope,:))
     180  call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
    171181ENDDO
    172182
    … …  
    176186
    177187  else !nslope=1 no slope, we copy all the values
    178     co2_ice_slope(:,:,1,:)=co2_ice_s(:,:,:)
    179     h2o_ice_s_slope(:,:,1,:)=h2o_ice_s(:,:,:)
    180     call get_var3("tsurf",tsurf_gcm(:,:,1,:))
     188
     189    CALL get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
     190    CALL get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
     191    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
     192#ifndef CPP_STD
     193    call get_var3("watercap", watercap_slope(:,:,1,:))
     194#endif
     195
    181196    if(soil_pem) then
    182       call get_var4("tsoil",tsoil_gcm(:,:,:,1,:))
     197      call get_var4("tsoil",tsoil_gcm_dyn(:,:,:,1,:))
    183198      call get_var4("inertiesoil",TI_gcm(:,:,:,1,:))
    184199    endif !soil_pem
    … …  
    187202! Compute the minimum over the year for each point
    188203  print *, "Computing the min of h2o_ice_slope"
    189 !  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope+watercap_slope,4)
    190   min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
     204  min_h2o_ice_dyn(:,:,:)=minval(h2o_ice_s_dyn+watercap_slope,4)
     205!  min_h2o_ice_dyn(:,:,:)=minval(h2o_ice_s_dyn,4)
    191206  print *, "Computing the min of co2_ice_slope"
    192   min_co2_ice_slope(:,:,:)=minval(co2_ice_slope,4)
     207  min_co2_ice_dyn(:,:,:)=minval(co2_ice_slope_dyn,4)
    193208
    194209!Compute averages
    195210
    196211    print *, "Computing average of tsurf"
    197     tsurf_ave(:,:,:)=SUM(tsurf_gcm(:,:,:,:),4)/timelen
     212    tsurf_ave_dyn(:,:,:)=SUM(tsurf_gcm_dyn(:,:,:,:),4)/timelen
     213
     214  DO islope = 1,nslope
     215    DO t=1,timelen
     216      CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
     217    ENDDO
     218  ENDDO
    198219
    199220  if(soil_pem) then
    200221    print *, "Computing average of tsoil"
    201     tsoil_ave(:,:,:,:)=SUM(tsoil_gcm(:,:,:,:,:),5)/timelen
     222    tsoil_ave_dyn(:,:,:,:)=SUM(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
    202223    print *, "Computing average of TI"
    203     TI_ave(:,:,:,:)=SUM(TI_gcm(:,:,:,:,:),5)/timelen
     224    TI_ave_dyn(:,:,:,:)=SUM(TI_gcm(:,:,:,:,:),5)/timelen
     225    print *, "Computing average of watersurf_density"
     226    watersurf_density_ave(:,:) = SUM(watersurf_density(:,:,:),3)/timelen
    204227  endif
    205228
    … …  
    208231    DO j = 1, jjm+1
    209232       DO islope=1,nslope
    210           if (min_co2_ice_slope(i,j,islope).LT.0) then
    211             min_co2_ice_slope(i,j,islope)  = 0.
     233          if (min_co2_ice_dyn(i,j,islope).LT.0) then
     234            min_co2_ice_dyn(i,j,islope)  = 0.
    212235          endif
    213           if (min_h2o_ice_slope(i,j,islope).LT.0) then
    214             min_h2o_ice_slope(i,j,islope)  = 0.
     236          if (min_h2o_ice_dyn(i,j,islope).LT.0) then
     237            min_h2o_ice_dyn(i,j,islope)  = 0.
    215238          endif
    216239       ENDDO
    … …  
    221244    DO j = 1, jjm+1
    222245      DO t = 1, timelen
    223          if (q_co2_GCM(i,j,t).LT.0) then
    224               q_co2_GCM(i,j,t)=1E-10
    225          elseif (q_co2_GCM(i,j,t).GT.1) then
    226               q_co2_GCM(i,j,t)=1.
     246         if (q_co2_dyn(i,j,t).LT.0) then
     247              q_co2_dyn(i,j,t)=1E-10
     248         elseif (q_co2_dyn(i,j,t).GT.1) then
     249              q_co2_dyn(i,j,t)=1.
    227250         endif
    228          if (q_h2o_GCM(i,j,t).LT.0) then
    229               q_h2o_GCM(i,j,t)=1E-30
    230          elseif (q_h2o_GCM(i,j,t).GT.1) then
    231               q_h2o_GCM(i,j,t)=1.
     251         if (q_h2o_dyn(i,j,t).LT.0) then
     252              q_h2o_dyn(i,j,t)=1E-30
     253         elseif (q_h2o_dyn(i,j,t).GT.1) then
     254              q_h2o_dyn(i,j,t)=1.
    232255         endif
    233          mmean=1/(A*q_co2_GCM(i,j,t) +B)
    234          vmr_co2_gcm(i,j,t) = q_co2_GCM(i,j,t)*mmean/m_co2
     256         mmean=1/(A*q_co2_dyn(i,j,t) +B)
     257         vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
    235258      ENDDO
    236259    ENDDO
    237260  ENDDO
    238261
    239       deallocate(h2o_ice_s_slope)
     262     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
     263     call gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,ps_GCM,ps_timeseries)
     264     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,q_co2_dyn,q_co2)
     265     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,q_h2o_dyn,q_h2o)
     266
     267     DO islope = 1,nslope
     268       CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
     269       CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
     270       if(soil_pem) then
     271         CALL gr_dyn_fi(nsoilmx,iim_input+1,jjm_input+1,ngrid,TI_ave_dyn(:,:,:,islope),TI_ave(:,:,islope))
     272       DO l=1,nsoilmx
     273         CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
     274         DO t=1,timelen
     275           CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
     276           CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
     277         ENDDO
     278       ENDDO
     279       endif !soil_pem
     280       DO t=1,timelen
     281         CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
     282         CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
     283       ENDDO
     284     ENDDO
     285
     286     CALL gr_dyn_fi(nslope,iim_input+1,jjm_input+1,ngrid,tsurf_ave_dyn,tsurf_ave)
    240287
    241288  CONTAINS
Note: See TracChangeset for help on using the changeset viewer.