Changeset 3076


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

PEM:
Big cleaning/improvements of the PEM:

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

JBC

Location:
trunk
Files:
12 edited
11 moved

Legend:

Unmodified
Added
Removed

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

    r3075 r3076  
    1 !
    2 ! $Id: abort_pem.F
    3 !
    4 c
    5 c
    6       SUBROUTINE abort_pem(modname, message, ierr)
    7      
     1MODULE abort_pem_mod
     2
     3implicit none
     4
     5!=======================================================================
     6contains
     7!=======================================================================
     8
     9SUBROUTINE abort_pem(modname,message,ierr)
     10
    811#ifdef CPP_IOIPSL
    9       USE IOIPSL
     12    use IOIPSL
    1013#else
    11 ! if not using IOIPSL, we still need to use (a local version of) getin_dump
    12       USE ioipsl_getincom
     14    ! if not using IOIPSL, we still need to use (a local version of) getin_dump
     15    use ioipsl_getincom
    1316#endif
    1417
    1518#ifdef CPP_XIOS
    16     ! ug Pour les sorties XIOS
    17       USE wxios
     19    use wxios ! ug Pour les sorties XIOS
    1820#endif
    1921
     22implicit none
     23
    2024#include "iniprint.h"
    21  
    22 C
    23 C Stops the simulation cleanly, closing files and printing various
    24 C comments
    25 C
    26 C  Input: modname = name of calling program
    27 C         message = stuff to print
    28 C         ierr    = severity of situation ( = 0 normal )
    2925
    30       character(len=*), intent(in):: modname
    31       integer, intent(in):: ierr
    32       character(len=*), intent(in):: message
     26! Stop the simulation cleanly, closing files and printing various comments
     27! Input: modname = name of calling program
     28!        message = stuff to print
     29!        ierr    = severity of situation (= 0 normal)
     30character(len = *), intent(in) :: modname
     31integer,            intent(in) :: ierr
     32character(len = *), intent(in) :: message
    3333
    34       write(lunout,*) 'in abort_pem'
     34!----- Code
     35write(lunout,*) 'in abort_pem'
    3536
    3637#ifdef CPP_XIOS
    37     !Fermeture propre de XIOS
    38       CALL wxios_close()
     38    CALL wxios_close() ! Fermeture propre de XIOS
    3939#endif
    4040
    4141#ifdef CPP_IOIPSL
    42       call histclo
    43       call restclo
     42    call histclo
     43    call restclo
    4444#endif
    45       call getin_dump
    46       write(lunout,*) 'Stopping in ', modname
    47       write(lunout,*) 'Reason = ',message
    48       if (ierr .eq. 0) then
    49         write(lunout,*) 'Everything is cool'
    50         stop
    51       else
    52         write(lunout,*) 'Houston, we have a problem ', ierr
    53         stop 1
    54       endif
    55       END
     45
     46call getin_dump
     47write(lunout,*) 'Stopping in ', modname
     48write(lunout,*) 'Reason = ', message
     49if (ierr == 0) then
     50    write(lunout,*) 'Everything is cool'
     51    stop
     52else
     53    write(lunout,*) 'Houston, we have a problem ', ierr
     54    stop 1
     55endif
     56
     57END SUBROUTINE abort_pem
     58
     59END MODULE abort_pem_mod

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3070 r3076  
    8989== 05/10/2023 == JBC
    9090Correction of a bug: the variable 'g' was not correctly initialized in 1D. A little of code cleaning.
     91
     92== 06/10/2023 == JBC
     93Big cleaning/improvements of the PEM:
     94    - Conversion of "abort_pem.F" and "soil_settings_PEM.F" into Fortran 90;
     95    - Transformation of every PEM subroutines into module;
     96    - Rewriting of many subroutines with modern Fortran syntax;
     97    - Correction of a bug in "pem.F90" when calling 'recomp_tend_co2_slope'. The arguments were given in disorder and emissivity was missing;
     98    - Update of "launch_pem.sh" in deftank.

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

    r3075 r3076  
     1MODULE compute_tendencies_slope_mod
     2
     3implicit none
     4
     5!=======================================================================
     6contains
     7!=======================================================================
     8
    19SUBROUTINE compute_tendencies_slope(ngrid,nslope,min_ice_Y1,min_ice_Y2,tendencies_ice)
    210
    … …  
    1826!   OUTPUT
    1927real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of perenial ice
    20 
    2128!=======================================================================
    2229
    … …  
    2936END SUBROUTINE compute_tendencies_slope
    3037
     38END MODULE compute_tendencies_slope_mod
     39

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

    r3039 r3076  
    88!=======================================================================
    99
    10 IMPLICIT NONE
     10implicit none
    1111
    12 CONTAINS
     12!=======================================================================
     13contains
     14!=======================================================================
    1315
    14   SUBROUTINE conf_pem(i_myear,n_myears)
     16SUBROUTINE conf_pem(i_myear,n_myears)
    1517
    1618#ifdef CPP_IOIPSL
    17   use IOIPSL, only: getin
     19    use IOIPSL,          only: getin
    1820#else
    19   ! if not using IOIPSL, we still need to use (a local version of) getin
    20   use ioipsl_getincom, only: getin
     21    ! if not using IOIPSL, we still need to use (a local version of) getin
     22    use ioipsl_getincom, only: getin
    2123#endif
    22  
    23   use time_evol_mod,  only: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
    24                             evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
    25   use comsoil_h_pem,  only: soil_pem, fluxgeo, water_reservoir_nom, depth_breccia, depth_bedrock, reg_thprop_dependp
    26   use adsorption_mod, only: adsorption_pem
    27   use glaciers_mod,   only: co2glaciersflow, h2oglaciersflow
    28   use ice_table_mod,  only: icetable_equilibrium, icetable_dynamic
    2924
    30   integer, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
     25use time_evol_mod,  only: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
     26                          evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
     27use comsoil_h_pem,  only: soil_pem, fluxgeo, water_reservoir_nom, depth_breccia, depth_bedrock, reg_thprop_dependp
     28use adsorption_mod, only: adsorption_pem
     29use glaciers_mod,   only: co2glaciersflow, h2oglaciersflow
     30use ice_table_mod,  only: icetable_equilibrium, icetable_dynamic
    3131
    32   character(len=20), parameter :: modname ='conf_pem'
    33   integer                      :: ierr
    34   integer                      :: year_earth_bp_ini ! Initial year (in Earth years) of the simulation of the PEM defined in run.def
     32implicit none
     33
     34integer, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
     35
     36character(20), parameter :: modname ='conf_pem'
     37integer                  :: ierr
     38integer                  :: year_earth_bp_ini ! Initial year (in Earth years) of the simulation of the PEM defined in run.def
    3539
    3640!PEM parameters
    … …  
    4852
    4953! #---------- Orbital parameters
    50   evol_orbit_pem=.false.
    51   call getin('evol_orbit_pem',evol_orbit_pem)
     54evol_orbit_pem = .false.
     55call getin('evol_orbit_pem',evol_orbit_pem)
    5256
    53   year_bp_ini = 0.
    54   call getin('year_earth_bp_ini',year_earth_bp_ini)
    55   year_bp_ini = year_earth_bp_ini/convert_years
     57year_bp_ini = 0.
     58call getin('year_earth_bp_ini',year_earth_bp_ini)
     59year_bp_ini = year_earth_bp_ini/convert_years
    5660
    57   var_obl = .true.
    58   call getin('var_obl',var_obl)
    59   write(*,*) 'Does obliquity vary ?',var_obl
     61var_obl = .true.
     62call getin('var_obl',var_obl)
     63write(*,*) 'Does obliquity vary ?',var_obl
    6064
    61   var_ecc = .true.
    62   call getin('var_ecc',var_ecc)
    63   write(*,*) 'Does excentricity vary ?',var_ecc
     65var_ecc = .true.
     66call getin('var_ecc',var_ecc)
     67write(*,*) 'Does excentricity vary ?',var_ecc
    6468
    65   var_lsp = .true.
    66   call getin('var_lsp',var_lsp)
    67   write(*,*) 'Does Ls peri vary ?',var_lsp
     69var_lsp = .true.
     70call getin('var_lsp',var_lsp)
     71write(*,*) 'Does Ls peri vary ?',var_lsp
    6872
    6973! #---------- Stopping criteria parameters
    70   Max_iter_pem=100000000
    71   call getin('Max_iter_pem', Max_iter_pem)
     74Max_iter_pem = 100000000
     75call getin('Max_iter_pem',Max_iter_pem)
    7276
    73   water_ice_criterion=0.2
    74   call getin('water_ice_criterion', water_ice_criterion)
     77water_ice_criterion = 0.2
     78call getin('water_ice_criterion',water_ice_criterion)
    7579
    76   co2_ice_criterion=0.2
    77   call getin('co2_ice_criterion', co2_ice_criterion)
     80co2_ice_criterion = 0.2
     81call getin('co2_ice_criterion',co2_ice_criterion)
    7882
    79   ps_criterion = 0.15
    80   call getin('ps_criterion',ps_criterion)
     83ps_criterion = 0.15
     84call getin('ps_criterion',ps_criterion)
    8185
    82   dt_pem=1
    83   call getin('dt_pem', dt_pem)
     86dt_pem = 1
     87call getin('dt_pem', dt_pem)
    8488
    8589! #---------- Subsurface parameters
    86   soil_pem=.true.
    87   call getin('soil_pem', soil_pem)
     90soil_pem = .true.
     91call getin('soil_pem',soil_pem)
    8892
    89   adsorption_pem = .true.
    90   call getin('adsorption_pem',adsorption_pem)
     93adsorption_pem = .true.
     94call getin('adsorption_pem',adsorption_pem)
    9195
    92   co2glaciersflow = .true.
    93   call getin('co2glaciersflow', co2glaciersflow)
     96co2glaciersflow = .true.
     97call getin('co2glaciersflow',co2glaciersflow)
    9498
    95   h2oglaciersflow = .true.
    96   call getin('h2oglaciersflow', h2oglaciersflow)
     99h2oglaciersflow = .true.
     100call getin('h2oglaciersflow',h2oglaciersflow)
    97101
    98   reg_thprop_dependp = .false.
    99   call getin('reg_thprop_dependp',reg_thprop_dependp)
    100   write(*,*)  'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
     102reg_thprop_dependp = .false.
     103call getin('reg_thprop_dependp',reg_thprop_dependp)
     104write(*,*)  'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
    101105
    102106! #---------- Layering parameters
    103   fluxgeo = 0.
    104   call getin('fluxgeo',fluxgeo)
    105   write(*,*) 'Flux Geothermal is set to',fluxgeo
    106    
    107   depth_breccia   = 10.
    108   call getin('depth_breccia',depth_breccia)
    109   write(*,*) 'Depth of breccia is set to',depth_breccia
     107fluxgeo = 0.
     108call getin('fluxgeo',fluxgeo)
     109write(*,*) 'Flux Geothermal is set to',fluxgeo
    110110
    111   depth_bedrock   = 1000.
    112   call getin('depth_bedrock',depth_bedrock)
    113   write(*,*) 'Depth of bedrock is set to',depth_bedrock
     111depth_breccia = 10.
     112call getin('depth_breccia',depth_breccia)
     113write(*,*) 'Depth of breccia is set to',depth_breccia
    114114
    115    icetable_equilibrium = .true.
    116    call getin('icetable_equilibrium',icetable_equilibrium)
    117    write(*,*)  'Is the ice table computed at equilibrium?', icetable_equilibrium
     115depth_bedrock = 1000.
     116call getin('depth_bedrock',depth_bedrock)
     117write(*,*) 'Depth of bedrock is set to',depth_bedrock
    118118
    119    icetable_dynamic = .false.
    120    call getin('icetable_dynamic',icetable_dynamic)
    121    write(*,*)  'Is the ice table computed with the dynamic method?', icetable_dynamic
    122   if ((.not.soil_pem).and.((icetable_equilibrium).or.(icetable_dynamic))) then
    123        write(*,*) 'Ice table  must be used when soil_pem = T'
    124        call abort_physic(modname,"Ice table  must be used when soil_pem = T",1)
    125   endif
     119icetable_equilibrium = .true.
     120call getin('icetable_equilibrium',icetable_equilibrium)
     121write(*,*)  'Is the ice table computed at equilibrium?', icetable_equilibrium
    126122
    127   if ((.not.soil_pem).and.adsorption_pem) then
    128        write(*,*) 'Adsorption must be used when soil_pem = T'
    129        call abort_physic(modname,"Adsorption must be used when soil_pem = T",1)
    130   endif
    131  
    132   if ((.not.soil_pem).and.(fluxgeo.gt.0.)) then
    133        write(*,*) 'Soil is not activated but Flux Geo > 0.'
    134        call abort_physic(modname,"Soil is not activated but Flux Geo > 0.",1)
    135   endif
    136  
    137   if ((.not.soil_pem).and.reg_thprop_dependp) then
    138      write(*,*) 'Regolith properties vary with Ps only when soil is set to true'
    139      call abort_physic(modname,'Regolith properties vary with Ps only when soil is set to true',1)
    140   endif
     123icetable_dynamic = .false.
     124call getin('icetable_dynamic',icetable_dynamic)
     125write(*,*)  'Is the ice table computed with the dynamic method?', icetable_dynamic
     126if ((.not. soil_pem) .and. ((icetable_equilibrium) .or. (icetable_dynamic))) then
     127    write(*,*) 'Ice table  must be used when soil_pem = T'
     128    call abort_physic(modname,"Ice table  must be used when soil_pem = T",1)
     129endif
    141130
    142   if (evol_orbit_pem .and. year_bp_ini == 0.) then
    143      write(*,*)  'You want to follow the file obl_ecc_lsp.asc for changing orb parameters,'
    144      write(*,*)  'but you did not specify from which year to start.'
    145      call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini=0",1)
    146   endif
     131if ((.not. soil_pem) .and. adsorption_pem) then
     132    write(*,*) 'Adsorption must be used when soil_pem = T'
     133    call abort_physic(modname,"Adsorption must be used when soil_pem = T",1)
     134endif
    147135
    148   water_reservoir_nom = 1e4
    149   call getin('water_reservoir_nom',water_reservoir_nom)
     136if ((.not. soil_pem) .and. (fluxgeo > 0.)) then
     137    write(*,*) 'Soil is not activated but Flux Geo > 0.'
     138    call abort_physic(modname,"Soil is not activated but Flux Geo > 0.",1)
     139endif
    150140
    151   END SUBROUTINE conf_pem
     141if ((.not. soil_pem) .and. reg_thprop_dependp) then
     142    write(*,*) 'Regolith properties vary with Ps only when soil is set to true'
     143    call abort_physic(modname,'Regolith properties vary with Ps only when soil is set to true',1)
     144endif
     145
     146if (evol_orbit_pem .and. year_bp_ini == 0.) then
     147    write(*,*)  'You want to follow the file obl_ecc_lsp.asc for changing orb parameters,'
     148    write(*,*)  'but you did not specify from which year to start.'
     149    call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini=0",1)
     150endif
     151
     152water_reservoir_nom = 1.e4
     153call getin('water_reservoir_nom',water_reservoir_nom)
     154
     155END SUBROUTINE conf_pem
    152156
    153157END MODULE conf_pem_mod

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

    r2998 r3076  
    1     MODULE constants_marspem_mod
     1MODULE constants_marspem_mod
    22
     3implicit none
    34
    4     IMPLICIT NONE
    55! Duration of a year and day
    6       INTEGER,PARAMETER :: sols_per_my =669 ! Number of Sols per year
    7       REAL, PARAMETER :: sec_per_sol=88775. ! Duration of a sol, in seconds
     6integer, parameter :: sols_per_my = 669    ! Number of Sols per year
     7real,    parameter :: sec_per_sol = 88775. ! Duration of a sol, in seconds
    88
    99! Molecular masses for CO2,H2O and non condensible gaz, following Franz et al. 2017   
    10       REAL,PARAMETER :: m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
    11       REAL,PARAMETER :: m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
    12       REAL,PARAMETER :: m_h2o = 18.01528E-3      ! Molecular weight of h2o (kg/mol)
     10real, parameter :: m_co2 = 44.01e-3    ! CO2 molecular mass (kg/mol)   
     11real, parameter :: m_noco2 = 33.37e-3  ! Non condensible mol mass (kg/mol)   
     12real, parameter :: m_h2o = 18.01528e-3 ! Molecular weight of h2o (kg/mol)
    1313
    1414!     Coefficient for Clapeyron law for CO2 condensation temperature (Tco2 = beta/(alpha-log(vmr)),following James et al. 1992
    15       REAL,PARAMETER :: alpha_clap_co2 = 23.3494  !Uniteless,James et al. 1992
    16       REAL,PARAMETER :: beta_clap_co2 = 3182.48   !Kelvin, James et al. 1992
     15real, parameter :: alpha_clap_co2 = 23.3494 ! Uniteless, James et al. 1992
     16real, parameter :: beta_clap_co2 = 3182.48  ! Kelvin, James et al. 1992
    1717
    1818!     Coefficient for Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
    19       REAL,PARAMETER :: alpha_clap_h2o = 28.9074 ! Uniteless, Murphy and Koop 2005
    20       REAL,PARAMETER :: beta_clap_h2o = -6143.7   ! Kelvin, Murphy and Koop 2005
     19real, parameter :: alpha_clap_h2o = 28.9074 ! Uniteless, Murphy and Koop 2005
     20real, parameter :: beta_clap_h2o = -6143.7  ! Kelvin, Murphy and Koop 2005
    2121
    2222!     Density of the regolith (Zent et al., 1995, Buhler and Piqueux 2021)     
    23       REAL,PARAMETER ::  rho_regolith = 2000.    ! kg/m^3
     23real, parameter :: rho_regolith = 2000. ! kg/m^3
    2424
    2525!     Average  Thermal inertia of the surface, breccia, bedrock, following Mellon et al., 2000., Wood et al., 2008
    26       REAL,PARAMETER :: TI_regolith_avg = 250.        ! Averaged of the observed thermal inertia for regolith following Mellon et al., 2000[SI]
    27       REAL,PARAMETER :: TI_breccia = 750.             ! Thermal inertia of Breccia following Wood 2009 [SI]
    28       REAL,PARAMETER :: TI_bedrock = 2300.            ! Thermal inertia of Bedrock following Wood 2009 [SI]
     26real, parameter :: TI_regolith_avg = 250. ! Averaged of the observed thermal inertia for regolith following Mellon et al., 2000[SI]
     27real, parameter :: TI_breccia = 750.      ! Thermal inertia of Breccia following Wood 2009 [SI]
     28real, parameter :: TI_bedrock = 2300.     ! Thermal inertia of Bedrock following Wood 2009 [SI]
    2929
    3030!     Porosity of the soil
    31       REAL,PARAMETER :: porosity = 0.4                ! porosity of the martian soil, correspond to the value for a random loose packing of monodiperse sphere (Scott, 1960)
     31real, parameter :: porosity = 0.4 ! porosity of the martian soil, correspond to the value for a random loose packing of monodiperse sphere (Scott, 1960)
    3232
    3333!     Stefan Boltzmann constant
    34       REAL,PARAMETER :: sigmaB=5.678E-8
     34real, parameter :: sigmaB = 5.678e-8
    3535
    3636!     Latent heat of CO2
    37       REAL,PARAMETER ::  Lco2 =  5.71E5                ! Pilorget and Forget 2016
     37real, parameter :: Lco2 =  5.71e5 ! Pilorget and Forget 2016
    3838
    3939! Conversion H2O/CO2 frost to perenial frost and vice versa
    40       REAL,PARAMETER :: threshold_water_frost2perenial = 1000. !~ 1m
    41       REAL,PARAMETER :: threshold_co2_frost2perenial = 3200.   !~ 2m
     40real, parameter :: threshold_water_frost2perenial = 1000. !~ 1m
     41real, parameter :: threshold_co2_frost2perenial = 3200.   !~ 2m
    4242
     43END MODULE constants_marspem_mod
    4344
    44      
    45     END MODULE constants_marspem_mod
    46 

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

    r3002 r3076  
    2626!!!
    2727!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    28 
    2928
    3029IMPLICIT NONE
    … …  
    105104subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
    106105
    107      USE comconst_mod, ONLY: pi,g
     106use comconst_mod,  only: pi,g
     107use abort_pem_mod, only: abort_pem
    108108
    109109!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    181181!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    182182
    183       USE comconst_mod, ONLY: pi
     183use comconst_mod,  only: pi
     184use abort_pem_mod, only: abort_pem
    184185
    185186

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

    r3075 r3076  
     1MODULE info_run_PEM_mod
     2
     3implicit none
     4
     5!=======================================================================
     6contains
     7!=======================================================================
     8
    19SUBROUTINE info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
    210
    … …  
    917!=======================================================================
    1018
    11   use time_evol_mod, only: convert_years
     19use time_evol_mod, only: convert_years
    1220
    13   IMPLICIT NONE
     21implicit none
    1422
    15   integer, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
    16   integer, intent(in) :: i_myear, n_myear ! Current simulated Martian year and maximum number of Martian years to be simulated
     23!----- Arguments
     24integer, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
     25integer, intent(in) :: i_myear, n_myear          ! Current simulated Martian year and maximum number of Martian years to be simulated
    1726
    18   logical :: ok
     27!----- Local variables
     28logical :: ok
    1929
     30!----- Code
    2031inquire(file = 'info_run_PEM.txt', exist = ok)
    2132if (ok) then
    … …  
    3243
    3344END SUBROUTINE info_run_PEM
     45
     46END MODULE info_run_PEM_mod

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

    r3075 r3076  
    1    subroutine interpolate_TIPEM_TIGCM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_PEM,TI_GCM)
     1MODULE interpolate_TIPEM_TIGCM_mod
     2
     3implicit none
     4
     5!=======================================================================
     6contains
     7!=======================================================================
     8
     9SUBROUTINE interpolate_TIPEM_TIGCM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_PEM,TI_GCM)
    210
    311!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    412!!!
    513!!! Purpose: Transfer the thermal inertia from the PEM vertical  grid to the GCM vertical grid
    6 !!! 
    7 !!! 
     14!!!
     15!!!
    816!!! Author: LL
    917!!!
    1018!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1119
    12 
    13       implicit none
    14 
     20implicit none
    1521
    1622!======================================================================
    … …  
    1824!  ---------
    1925!  inputs:
    20       integer,intent(in) :: ngrid                                ! # of horizontal grid points
    21       integer,intent(in) :: nslope                               ! # of subslope wihtin the mesh
    22       integer,intent(in) :: nsoil_PEM                            ! # of soil layers in the PEM
    23       integer,intent(in) :: nsoil_GCM                            ! # of soil layers in the GCM
    24       real,intent(in) :: TI_PEM(ngrid,nsoil_PEM,nslope)          !  Thermal inertia in the PEM vertical grid [J/m^2/K/s^{1/2}]
    25       real,intent(inout) :: TI_GCM(ngrid,nsoil_GCM,nslope)       ! Thermal inertia in the GCM vertical grid [J/m^2/K/s^{1/2}]
     26integer,                                 intent(in) :: ngrid     ! # of horizontal grid points
     27integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
     28integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
     29integer,                                 intent(in) :: nsoil_GCM ! # of soil layers in the GCM
     30real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM    ! Thermal inertia in the PEM vertical grid [J/m^2/K/s^{1/2}]
    2631
    27 !local variable
    28       integer :: ig,islope,iloop                                 ! loop variables
    29      do ig = 1,ngrid
    30        do islope = 1,nslope
    31          do iloop = 1,nsoil_GCM
    32            TI_GCM(ig,iloop,islope) = TI_PEM(ig,iloop,islope)
    33          enddo
    34        enddo
    35      enddo
     32real, dimension(ngrid,nsoil_GCM,nslope), intent(inout) :: TI_GCM ! Thermal inertia in the GCM vertical grid [J/m^2/K/s^{1/2}]
    3633
     34!----- Code
     35TI_GCM(:,:,:) = TI_PEM(:,:nsoil_GCM,:)
    3736
    38   end subroutine interpolate_TIPEM_TIGCM
     37END SUBROUTINE interpolate_TIPEM_TIGCM
     38
     39END MODULE interpolate_TIPEM_TIGCM_mod

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

    r3039 r3076  
    1 module lask_param_mod
     1MODULE lask_param_mod
     2
    23!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    34!!!
    45!!! Purpose: Define parameters from Laskar et al., 2004 evolution
    5 !!! 
     6!!!
    67!!! Author: RV, JBC
    78!!!
    … …  
    1011implicit none
    1112
    12   real,save,allocatable :: yearlask(:) ! year before present from Laskar et al. Tab
    13   real,save,allocatable :: obllask(:)  ! obliquity    [deg]
    14   real,save,allocatable :: ecclask(:)  ! excentricity [deg]
    15   real,save,allocatable :: lsplask(:)  ! ls perihelie [deg]
    16   integer, save         :: last_ilask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
     13real, dimension(:), allocatable, save :: yearlask  ! year before present from Laskar et al. Tab
     14real, dimension(:), allocatable, save :: obllask    ! obliquity    [deg]
     15real, dimension(:), allocatable, save :: ecclask    ! excentricity [deg]
     16real, dimension(:), allocatable, save :: lsplask    ! ls perihelie [deg]
     17integer,                         save :: last_ilask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
    1718
     19!=======================================================================
    1820contains
     21!=======================================================================
    1922
    20   subroutine ini_lask_param_mod
    21  
    22   implicit none
     23SUBROUTINE ini_lask_param_mod
    2324
    24   integer :: nlask ! number of lines in Laskar files
    25   integer :: ierr
     25implicit none
    2626
    27   nlask = 0
    28   open(1,file = 'obl_ecc_lsp.asc')
    29   do
    30       read(1,*,iostat = ierr)
    31       if (ierr /= 0) exit
    32       nlask = nlask + 1
    33   enddo
    34   close(1)
    35   allocate(yearlask(nlask))
    36   allocate(obllask(nlask))
    37   allocate(ecclask(nlask))
    38   allocate(lsplask(nlask))
    39  
    40   end subroutine ini_lask_param_mod
     27integer :: nlask ! number of lines in Laskar files
     28integer :: ierr
    4129
     30nlask = 0
     31open(1,file = 'obl_ecc_lsp.asc')
     32do
     33    read(1,*,iostat = ierr)
     34    if (ierr /= 0) exit
     35    nlask = nlask + 1
     36enddo
     37close(1)
     38allocate(yearlask(nlask),obllask(nlask),ecclask(nlask),lsplask(nlask))
    4239
    43   subroutine end_lask_param_mod
     40END SUBROUTINE ini_lask_param_mod
    4441
    45   implicit none
     42!=======================================================================
    4643
    47   if (allocated(yearlask)) deallocate(yearlask)
    48   if (allocated(obllask)) deallocate(obllask)
    49   if (allocated(ecclask)) deallocate(ecclask)
    50   if (allocated(lsplask)) deallocate(lsplask)
     44SUBROUTINE end_lask_param_mod
    5145
    52   end subroutine end_lask_param_mod
    53  
    54 end module lask_param_mod
     46implicit none
     47
     48if (allocated(yearlask)) deallocate(yearlask)
     49if (allocated(obllask)) deallocate(obllask)
     50if (allocated(ecclask)) deallocate(ecclask)
     51if (allocated(lsplask)) deallocate(lsplask)
     52
     53END SUBROUTINE end_lask_param_mod
     54
     55END MODULE lask_param_mod

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

    r3075 r3076  
    1 !
    2 ! $Id $
    3 !
     1MODULE nb_time_step_GCM_mod
     2
     3use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
     4                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,       &
     5                  nf90_inquire_dimension, nf90_close
     6
     7implicit none
     8
     9character(256) :: msg, var, modname
     10
     11!=======================================================================
     12contains
     13!=======================================================================
     14
    415SUBROUTINE nb_time_step_GCM(fichnom,timelen)
    516
    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
     17implicit none
    1118
    1219!=======================================================================
    … …  
    1724!=======================================================================
    1825
    19   include "dimensions.h"
     26include "dimensions.h"
    2027
    21 !===============================================================================
     28!=======================================================================
    2229! Arguments:
    23   CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
    24 !===============================================================================
    25 !   Local Variables
    26   CHARACTER(LEN=256) :: msg, var, modname
    27   INTEGER :: iq, fID, vID, idecal
    28   INTEGER :: ierr
    29 
    30   INTEGER :: timelen ! number of times stored in the file
     30character(len = *), intent(in) :: fichnom !--- FILE NAME
     31!=======================================================================
     32!   Local Variables
     33integer        :: iq, fID, vID, idecal, ierr
     34integer        :: timelen ! number of times stored in the file
    3135!-----------------------------------------------------------------------
    32   modname="nb_time_step_GCM"
     36modname = "nb_time_step_GCM"
    3337
    3438!  Open initial state NetCDF file
    35   var=fichnom
    36   CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
     39var = fichnom
     40call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
    3741
    38       ierr = nf90_inq_varid (fID, "temps", vID)
    39       IF (ierr .NE. nf90_noerr) THEN
    40         write(*,*)"read_data_GCM: Le champ <temps> est absent"
    41         write(*,*)"read_data_GCM: J essaie <time_counter>"
    42         ierr = nf90_inq_varid (fID, "time_counter", vID)
    43         IF (ierr .NE. nf90_noerr) THEN
    44           write(*,*)"read_data_GCM: Le champ <time_counter> est absent"
    45           write(*,*)"read_data_GCM: J essaie <Time>"
    46           ierr = nf90_inq_varid (fID, "Time", vID)
    47           IF (ierr .NE. nf90_noerr) THEN
     42ierr = nf90_inq_varid (fID, "temps", vID)
     43if (ierr /= nf90_noerr) then
     44    write(*,*)"read_data_GCM: Le champ <temps> est absent"
     45    write(*,*)"read_data_GCM: J essaie <time_counter>"
     46    ierr = nf90_inq_varid (fID, "time_counter", vID)
     47    if (ierr /= nf90_noerr) then
     48        write(*,*)"read_data_GCM: Le champ <time_counter> est absent"
     49        write(*,*)"read_data_GCM: J essaie <Time>"
     50        ierr = nf90_inq_varid (fID, "Time", vID)
     51        if (ierr /= nf90_noerr) then
    4852            write(*,*)"read_data_GCM: Le champ <Time> est absent"
    4953            write(*,*)trim(nf90_strerror(ierr))
    50             CALL ABORT_gcm("nb_time_step_GCM", "", 1)
    51           ENDIF
    52           ! Get the length of the "Time" dimension
    53           ierr = nf90_inq_dimid(fID,"Time",vID)
    54           ierr = nf90_inquire_dimension(fID,vID,len=timelen)         
    55         ELSE
     54            call abort_gcm("nb_time_step_GCM", "", 1)
     55        endif
     56        ! Get the length of the "Time" dimension
     57        ierr = nf90_inq_dimid(fID,"Time",vID)
     58        ierr = nf90_inquire_dimension(fID,vID,len = timelen)
     59    else
    5660        ! Get the length of the "time_counter" dimension
    5761        ierr = nf90_inq_dimid(fID,"time_counter",vID)
    58         ierr = nf90_inquire_dimension(fID,vID,len=timelen)
    59         ENDIF
    60       ELSE   
    61         ! Get the length of the "temps" dimension
    62         ierr = nf90_inq_dimid(fID,"temps",vID)
    63         ierr = nf90_inquire_dimension(fID,vID,len=timelen)
    64       ENDIF
     62        ierr = nf90_inquire_dimension(fID,vID,len = timelen)
     63    endif
     64else
     65    ! Get the length of the "temps" dimension
     66    ierr = nf90_inq_dimid(fID,"temps",vID)
     67    ierr = nf90_inquire_dimension(fID,vID,len = timelen)
     68endif
    6569
    66   CALL err(NF90_CLOSE(fID),"close",fichnom)
     70call error_msg(NF90_CLOSE(fID),"close",fichnom)
    6771
    68   write(*,*) "The number of timestep of the PCM run data=", timelen
    69 
    70   CONTAINS
    71 
    72 SUBROUTINE err(ierr,typ,nam)
    73   INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
    74   CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
    75   CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
    76   IF(ierr==NF90_NoERR) RETURN
    77   SELECT CASE(typ)
    78     CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
    79     CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
    80     CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
    81     CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
    82   END SELECT
    83   CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
    84 END SUBROUTINE err
     72write(*,*) "The number of timestep of the PCM run data=", timelen
    8573
    8674END SUBROUTINE nb_time_step_GCM
     75
     76!=======================================================================
     77
     78SUBROUTINE error_msg(ierr,typ,nam)
     79
     80implicit none
     81
     82integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
     83character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
     84character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
     85
     86if (ierr == nf90_noerr) return
     87select case(typ)
     88    case('inq');   msg = "Field <"//trim(nam)//"> is missing"
     89    case('get');   msg = "Reading failed for <"//trim(nam)//">"
     90    case('open');  msg = "File opening failed for <"//trim(nam)//">"
     91    case('close'); msg = "File closing failed for <"//trim(nam)//">"
     92    case default
     93        write(*,*) 'There is no message for this error.'
     94        error stop
     95end select
     96call abort_gcm(trim(modname),trim(msg),ierr)
     97
     98END SUBROUTINE error_msg
     99
     100END MODULE nb_time_step_GCM_mod

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

    r3069 r3076  
    33!=======================================================================
    44!
    5 ! Purpose: Compute the maximum number of iteration for PEM based on the 
     5! Purpose: Compute the maximum number of iteration for PEM based on the
    66! stopping criterion given by the orbital parameters
    77!
    … …  
    99!=======================================================================
    1010
    11 IMPLICIT NONE
    12 
    13 CONTAINS
     11implicit none
     12
     13!=======================================================================
     14contains
     15!=======================================================================
    1416
    1517SUBROUTINE orbit_param_criterion(i_myear,year_iter_max)
    … …  
    2628    use planete_mod, only: e_elips, obliquit, lsperi
    2729#endif
    28   use time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
    29   use comconst_mod,   only: pi
    30   use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask
    31 
    32   IMPLICIT NONE
     30use time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
     31use comconst_mod,   only: pi
     32use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask
     33
     34implicit none
    3335
    3436!--------------------------------------------------------
    3537! Input Variables
    3638!--------------------------------------------------------
    37   integer, intent(in) :: i_myear ! Martian year of the simulation
     39integer, intent(in) :: i_myear ! Martian year of the simulation
    3840
    3941!--------------------------------------------------------
    4042! Output Variables
    4143!--------------------------------------------------------
    42   integer, intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much
    43 
    44 !--------------------------------------------------------
    45 ! Local variables 
    46 !--------------------------------------------------------
    47   integer :: Year                                           ! Year of the simulation
    48   real    :: Lsp                                            ! Ls perihelion
    49   real    :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
    50   real    :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
    51   real    :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
    52   integer :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
    53   integer :: ilask, iylask                                  ! Loop variables
    54   real    :: xa, xb, ya, yb                                 ! Variables for interpolation
    55   logical :: found_obl, found_ecc, found_lsp                ! Flag variables for orbital parameters
    56   logical :: crossed_modulo_d, crossed_modulo_i             ! Flag variables for modulo of Lsp
     44integer, intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much
     45
     46!--------------------------------------------------------
     47! Local variables
     48!--------------------------------------------------------
     49integer :: Year                                           ! Year of the simulation
     50real    :: Lsp                                            ! Ls perihelion
     51real    :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
     52real    :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
     53real    :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
     54integer :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
     55integer :: ilask, iylask                                  ! Loop variables
     56real    :: xa, xb, ya, yb                                 ! Variables for interpolation
     57logical :: found_obl, found_ecc, found_lsp                ! Flag variables for orbital parameters
     58logical :: crossed_modulo_d, crossed_modulo_i             ! Flag variables for modulo of Lsp
    5759
    5860! **********************************************************************
    … …  
    8486endif
    8587
    86 ! Constant max change case 
     88! Constant max change case
    8789max_change_obl = 0.6
    8890max_change_ecc = 2.8e-3
    … …  
    106108write(*,*) 'Maximum Lsp accepted =', max_lsp
    107109write(*,*) 'Minimum Lsp accepted =', min_lsp
    108 ! End Constant max change case 
     110! End Constant max change case
    109111
    110112! If we do not want some orb parameter to change, they should not be a stopping criterion,
    … …  
    175177            endif
    176178        else
    177             if (crossed_modulo_d) then ! If decreasing Lsp "crossed" the modulo but min_lsp has not been met yet 
     179            if (crossed_modulo_d) then ! If decreasing Lsp "crossed" the modulo but min_lsp has not been met yet
    178180                ya = ya - 360.
    179181                yb = yb - 360.
    … …  
    208210END SUBROUTINE orbit_param_criterion
    209211
    210 !********************************************************************************   
     212!***********************************************************************
    211213
    212214END MODULE orbit_param_criterion_mod

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

    r3072 r3076  
    2727PROGRAM pem
    2828
    29 use phyetat0_mod,               only: phyetat0
    30 use phyredem,                   only: physdem0, physdem1
    31 use netcdf,                     only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
    32 use turb_mod,                   only: q2, wstar
    33 use comslope_mod,               only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
    34 use logic_mod,                  only: iflag_phys
    35 use mod_const_mpi,              only: COMM_LMDZ
     29use phyetat0_mod,                 only: phyetat0
     30use phyredem,                     only: physdem0, physdem1
     31use netcdf,                       only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
     32use turb_mod,                     only: q2, wstar
     33use comslope_mod,                 only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
     34use logic_mod,                    only: iflag_phys
     35use mod_const_mpi,                only: COMM_LMDZ
    3636use infotrac
    37 use geometry_mod,               only: latitude_deg
    38 use conf_pem_mod,               only: conf_pem
    39 use pemredem,                   only: pemdem0, pemdem1
    40 use glaciers_mod,               only: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
    41 use criterion_pem_stop_mod,     only: criterion_waterice_stop, criterion_co2_stop
    42 use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
    43                                       m_noco2, threshold_water_frost2perenial, threshold_co2_frost2perenial
    44 use evol_co2_ice_s_mod,         only: evol_co2_ice_s
    45 use evol_h2o_ice_s_mod,         only: evol_h2o_ice_s
    46 use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
    47                                       TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
    48                                       tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
    49                                       fluxgeo,                          & ! Geothermal flux for the PEM and GCM
    50                                       water_reservoir                     ! Water ressources
    51 use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
    52                                       ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
    53                                       co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
    54 use time_evol_mod,              only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
    55 use orbit_param_criterion_mod,  only: orbit_param_criterion
    56 use recomp_orb_param_mod,       only: recomp_orb_param
    57 use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
    58                                       ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
    59 use soil_thermalproperties_mod, only: update_soil_thermalproperties
    60 use time_phylmdz_mod,           only: daysec, dtphys, day_end
     37use geometry_mod,                 only: latitude_deg
     38use conf_pem_mod,                 only: conf_pem
     39use pemredem,                     only: pemdem0, pemdem1
     40use glaciers_mod,                 only: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
     41use criterion_pem_stop_mod,       only: criterion_waterice_stop, criterion_co2_stop
     42use constants_marspem_mod,        only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
     43                                        m_noco2, threshold_water_frost2perenial, threshold_co2_frost2perenial
     44use evol_co2_ice_s_mod,           only: evol_co2_ice_s
     45use evol_h2o_ice_s_mod,           only: evol_h2o_ice_s
     46use comsoil_h_PEM,                only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
     47                                        TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
     48                                        tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
     49                                        fluxgeo,                          & ! Geothermal flux for the PEM and GCM
     50                                        water_reservoir                     ! Water ressources
     51use adsorption_mod,               only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
     52                                        ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
     53                                        co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
     54use time_evol_mod,                only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
     55use orbit_param_criterion_mod,    only: orbit_param_criterion
     56use recomp_orb_param_mod,         only: recomp_orb_param
     57use ice_table_mod,                only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
     58                                        ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
     59use soil_thermalproperties_mod,   only: update_soil_thermalproperties
     60use time_phylmdz_mod,             only: daysec, dtphys, day_end
     61use abort_pem_mod,                only: abort_pem
     62use soil_settings_PEM_mod,        only: soil_settings_PEM
     63use compute_tendencies_slope_mod, only: compute_tendencies_slope
     64use info_run_PEM_mod,             only: info_run_PEM
     65use interpolate_TIPEM_TIGCM_mod,  only: interpolate_TIPEM_TIGCM
     66use nb_time_step_GCM_mod,         only: nb_time_step_GCM
     67use pemetat0_mod,                 only: pemetat0
     68use read_data_GCM_mod,            only: read_data_GCM
     69use recomp_tend_co2_slope_mod,    only: recomp_tend_co2_slope
     70use soil_pem_compute_mod,         only: soil_pem_compute
    6171
    6272#ifndef CPP_STD
    … …  
    8595
    8696#ifndef CPP_1D
    87     use iniphysiq_mod,    only: iniphysiq
    88     use control_mod,      only: iphysiq, day_step, nsplit_phys
    89     use comconst_mod,     only: rad, g, cpp, pi
     97    use iniphysiq_mod,            only: iniphysiq
     98    use control_mod,              only: iphysiq, day_step, nsplit_phys
     99    use comconst_mod,             only: rad, g, cpp, pi
    90100#else
    91     use time_phylmdz_mod, only: iphysiq, day_step
    92     use comcstfi_h,       only: pi, rad, g, cpp
    93 #endif
    94 
    95 #ifdef CPP_1D
     101    use time_phylmdz_mod,         only: iphysiq, day_step
     102    use comcstfi_h,               only: pi, rad, g, cpp
    96103    use regular_lonlat_mod,       only: init_regular_lonlat
    97104    use physics_distribution_mod, only: init_physics_distribution
    … …  
    99106    use init_testphys1d_mod,      only: init_testphys1d
    100107    use comvert_mod,              only: ap, bp
     108    use writerestart1D_mod,       only: writerestart1D
    101109#endif
    102110
    103 IMPLICIT NONE
     111implicit none
    104112
    105113include "dimensions.h"
    … …  
    858866                Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
    859867                Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
    860                 call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
    861                 call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
     868                call soil_pem_compute(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
     869                call soil_pem_compute(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
    862870                tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
    863871                do ig = 1,ngrid
    … …  
    916924!------------------------
    917925    write(*,*) "Adaptation of the new co2 tendencies to the current pressure"
    918     call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
    919                                global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
     926    call recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),emis,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
     927                               global_ave_press_GCM,global_ave_press_new)
    920928
    921929!------------------------

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

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

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

    r3075 r3076  
    1 SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries, &
    2              min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
    3              watersurf_density_ave,watersoil_density)
    4 
    5       use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
    6                         nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
    7                         nf90_inquire_dimension,nf90_close
    8       use comsoil_h, only: nsoilmx
    9       USE comsoil_h_PEM, ONLY: soil_pem
    10       use constants_marspem_mod,only: m_co2,m_noco2
    11       IMPLICIT NONE
     1MODULE read_data_GCM_mod
     2
     3use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
     4                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,       &
     5                  nf90_inquire_dimension, nf90_close
     6
     7implicit none
     8
     9character(256) :: msg, var, modname ! for reading
     10integer        :: fID, vID          ! for reading
     11
     12!=======================================================================
     13contains
     14!=======================================================================
     15
     16SUBROUTINE read_data_GCM(fichnom,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,           &
     17                         min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
     18                         watersurf_density_ave,watersoil_density)
     19use comsoil_h,             only: nsoilmx
     20use comsoil_h_PEM,         only: soil_pem
     21use constants_marspem_mod, only: m_co2, m_noco2
     22
     23implicit none
    1224
    1325!=======================================================================
    … …  
    1830!=======================================================================
    1931
    20   include "dimensions.h"
    21 
     32include "dimensions.h"
     33
     34!=======================================================================
     35! Arguments:
     36
     37character(len = *), intent(in) :: fichnom                             !--- FILE NAME
     38integer,            intent(in) :: timelen                             ! number of times stored in the file
     39integer                        :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
     40! Ouputs
     41real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice  per slope of the year [kg/m^2]
     42real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
     43real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_gcm_phys ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
     44real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
     45real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
     46real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
     47real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per  slope of the year [kg/m^2]
     48!SOIL
     49real, dimension(ngrid,nslope),                 intent(out) :: tsurf_ave             ! Average surface temperature of the concatenated file [K]
     50real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_ave             ! Average soil temperature of the concatenated file [K]
     51real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_gcm             ! Surface temperature of the concatenated file, time series [K]
     52real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_gcm             ! Soil temperature of the concatenated file, time series [K]
     53real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! Water density at the surface [kg/m^3]
     54real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
    2255!===============================================================================
    23 ! Arguments:
    24 
    25   CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
    26   INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
    27   INTEGER :: iim_input,jjm_input,ngrid,nslope      ! number of points in the lat x lon dynamical grid, number of subgrid slopes
    28 ! Ouputs
    29   REAL, INTENT(OUT) ::  min_co2_ice(ngrid,nslope) ! Minimum of co2 ice  per slope of the year [kg/m^2]
    30   REAL, INTENT(OUT) ::  min_h2o_ice(ngrid,nslope) ! Minimum of h2o ice per slope of the year [kg/m^2]
    31   REAL, INTENT(OUT)  :: vmr_co2_gcm_phys(ngrid,timelen) ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
    32   REAL, INTENT(OUT) ::  ps_timeseries(ngrid,timelen)! Surface Pressure [Pa]
    33   REAL, INTENT(OUT) ::  q_co2(ngrid,timelen)        ! CO2 mass mixing ratio in the first layer [kg/m^3]
    34   REAL, INTENT(OUT) ::  q_h2o(ngrid,timelen)        ! H2O mass mixing ratio in the first layer [kg/m^3]
    35   REAL, INTENT(OUT) ::  co2_ice_slope(ngrid,nslope,timelen) ! co2 ice amount per  slope of the year [kg/m^2]
    36 !SOIL
    37   REAL, INTENT(OUT) ::  tsurf_ave(ngrid,nslope)         ! Average surface temperature of the concatenated file [K]
    38   REAL, INTENT(OUT) ::  tsoil_ave(ngrid,nsoilmx,nslope) ! Average soil temperature of the concatenated file [K]
    39   REAL ,INTENT(OUT) ::  tsurf_gcm(ngrid,nslope,timelen)                  ! Surface temperature of the concatenated file, time series [K]
    40   REAL , INTENT(OUT) ::  tsoil_gcm(ngrid,nsoilmx,nslope,timelen)         ! Soil temperature of the concatenated file, time series [K]
    41   REAL , INTENT(OUT) ::  watersurf_density_ave(ngrid,nslope)             ! Water density at the surface [kg/m^3]
    42   REAL , INTENT(OUT) ::  watersoil_density(ngrid,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
    43 !===============================================================================
    44 !   Local Variables
    45   CHARACTER(LEN=256) :: msg, var, modname               ! for reading
    46   INTEGER :: iq, fID, vID, idecal                       ! for reading
    47   INTEGER :: ierr                                       ! for reading
    48   CHARACTER(len=12) :: start_file_type="earth" ! default start file type
    49 
    50   REAL,ALLOCATABLE :: time(:) ! times stored in start
    51   INTEGER :: indextime ! index of selected time
    52 
    53   INTEGER :: edges(4),corner(4)
    54   INTEGER :: i,j,l,t                                                   ! loop variables
    55   real  ::  A , B, mmean                                           ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
    56 
    57   INTEGER :: islope                                                    ! loop for variables
    58   CHARACTER*2 :: num                                                   ! for reading sloped variables
    59   REAL ::  h2o_ice_s_dyn(iim_input+1,jjm_input+1,nslope,timelen)       ! h2o ice per slope of the concatenated file [kg/m^2]
    60   REAL ::  watercap_slope(iim_input+1,jjm_input+1,nslope,timelen)
    61   REAL ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)                ! CO2 volume mixing ratio in the first layer  [mol/m^3]
    62   REAL ::  ps_GCM(iim_input+1,jjm_input+1,timelen)                     ! Surface Pressure [Pa]
    63   REAL ::  min_co2_ice_dyn(iim_input+1,jjm_input+1,nslope)
    64   REAL ::  min_h2o_ice_dyn(iim_input+1,jjm_input+1,nslope)
    65   REAL ::  tsurf_ave_dyn(iim_input+1,jjm_input+1,nslope)               ! Average surface temperature of the concatenated file [K]
    66   REAL ::  tsoil_ave_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope)       ! Average soil temperature of the concatenated file [K]
    67   REAL ::  tsurf_gcm_dyn(iim_input+1,jjm_input+1,nslope,timelen)       ! Surface temperature of the concatenated file, time series [K]
    68   REAL ::  tsoil_gcm_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)! Soil temperature of the concatenated file, time series [K]
    69   REAL ::  q_co2_dyn(iim_input+1,jjm_input+1,timelen)                  ! CO2 mass mixing ratio in the first layer [kg/m^3]
    70   REAL ::  q_h2o_dyn(iim_input+1,jjm_input+1,timelen)                  ! H2O mass mixing ratio in the first layer [kg/m^3]
    71   REAL ::  co2_ice_slope_dyn(iim_input+1,jjm_input+1,nslope,timelen)  ! co2 ice amount per  slope of the year [kg/m^2]
    72   REAL ::  watersurf_density_dyn(iim_input+1,jjm_input+1,nslope,timelen)! Water density at the surface, time series [kg/m^3]
    73   REAL ::  watersurf_density(ngrid,nslope,timelen)                     ! Water density at the surface, time series [kg/m^3]
    74   REAL ::  watersoil_density_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
     56!   Local Variables
     57character(12)                   :: start_file_type = "earth" ! default start file type
     58real, dimension(:), allocatable :: time      ! times stored in start
     59integer                         :: indextime ! index of selected time
     60integer                         :: edges(4), corner(4)
     61integer                         :: i, j, l, t          ! loop variables
     62real                            ::  A , B, mmean       ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
     63integer                         :: islope              ! loop for variables
     64character(2)                    :: num                 ! for reading sloped variables
     65real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn         ! h2o ice per slope of the concatenated file [kg/m^2]
     66real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap_slope
     67real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm           ! CO2 volume mixing ratio in the first layer  [mol/m^3]
     68real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                ! Surface Pressure [Pa]
     69real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
     70real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
     71real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn         ! Average surface temperature of the concatenated file [K]
     72real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn         ! Average soil temperature of the concatenated file [K]
     73real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn         ! Surface temperature of the concatenated file, time series [K]
     74real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn         ! Soil temperature of the concatenated file, time series [K]
     75real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn             ! CO2 mass mixing ratio in the first layer [kg/m^3]
     76real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn             ! H2O mass mixing ratio in the first layer [kg/m^3]
     77real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn     ! co2 ice amount per  slope of the year [kg/m^2]
     78real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn ! Water density at the surface, time series [kg/m^3]
     79real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn ! Water density in the soil layer, time series [kg/m^3]
     80real, dimension(ngrid,nslope,timelen)                               :: watersurf_density     ! Water density at the surface, time series [kg/m^3]
    7581
    7682!-----------------------------------------------------------------------
    77   modname="read_data_gcm"
    78 
    79       A =(1/m_co2 - 1/m_noco2)
    80       B=1/m_noco2
    81 
    82   write(*,*) "Opening ", fichnom, "..."
     83modname="read_data_gcm"
     84
     85A = (1/m_co2 - 1/m_noco2)
     86B = 1/m_noco2
     87
     88write(*,*) "Opening ", fichnom, "..."
    8389
    8490!  Open initial state NetCDF file
    85   var=fichnom
    86   CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
    87 
    88      write(*,*) "Downloading data for vmr co2..."
    89 
    90   CALL get_var3("co2_layer1"   ,q_co2_dyn)
    91 
    92      write(*,*) "Downloading data for vmr co2 done"
    93      write(*,*) "Downloading data for vmr h20..."
    94 
    95   CALL get_var3("h2o_layer1"   ,q_h2o_dyn)
    96 
    97      write(*,*) "Downloading data for vmr h2o done"
    98      write(*,*) "Downloading data for surface pressure ..."
    99 
    100   CALL get_var3("ps"   ,ps_GCM)
    101 
    102      write(*,*) "Downloading data for surface pressure done"
    103      write(*,*) "nslope=", nslope
    104 
    105 if(nslope.gt.1) then
    106 
    107      write(*,*) "Downloading data for co2ice_slope ..."
    108 
    109 DO islope=1,nslope
    110   write(num,fmt='(i2.2)') islope
    111   call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
    112 ENDDO
    113 
    114      write(*,*) "Downloading data for co2ice_slope done"
    115      write(*,*) "Downloading data for h2o_ice_s_slope ..."
    116 
    117 DO islope=1,nslope
    118   write(num,fmt='(i2.2)') islope
    119   call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
    120 ENDDO
    121 
    122      write(*,*) "Downloading data for h2o_ice_s_slope done"
     91var = fichnom
     92call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
     93
     94write(*,*) "Downloading data for vmr co2..."
     95
     96call get_var3("co2_layer1"   ,q_co2_dyn)
     97
     98write(*,*) "Downloading data for vmr co2 done"
     99write(*,*) "Downloading data for vmr h20..."
     100
     101call get_var3("h2o_layer1"   ,q_h2o_dyn)
     102
     103write(*,*) "Downloading data for vmr h2o done"
     104write(*,*) "Downloading data for surface pressure ..."
     105
     106call get_var3("ps"   ,ps_GCM)
     107
     108write(*,*) "Downloading data for surface pressure done"
     109write(*,*) "nslope=", nslope
     110
     111if (nslope > 1) then
     112    write(*,*) "Downloading data for co2ice_slope ..."
     113
     114    do islope = 1,nslope
     115        write(num,'(i2.2)') islope
     116        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
     117    enddo
     118   
     119    write(*,*) "Downloading data for co2ice_slope done"
     120    write(*,*) "Downloading data for h2o_ice_s_slope ..."
     121   
     122    do islope = 1,nslope
     123        write(num,'(i2.2)') islope
     124        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
     125    enddo
     126
     127    write(*,*) "Downloading data for h2o_ice_s_slope done"
    123128
    124129#ifndef CPP_STD
    125 
    126      write(*,*) "Downloading data for watercap_slope ..."
    127 DO islope=1,nslope
    128        write(num,fmt='(i2.2)') islope
    129        call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
    130 !        watercap_slope(:,:,:,:)= 0.
    131 ENDDO           
    132      write(*,*) "Downloading data for watercap_slope done"
    133 
    134 #endif
     130    write(*,*) "Downloading data for watercap_slope ..."
     131    do islope = 1,nslope
     132        write(num,'(i2.2)') islope
     133        call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
     134!        watercap_slope(:,:,:,:) = 0.
     135    enddo
     136    write(*,*) "Downloading data for watercap_slope done"
     137#endif
     138
     139    write(*,*) "Downloading data for tsurf_slope ..."
     140
     141    do islope=1,nslope
     142        write(num,'(i2.2)') islope
     143        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
     144    enddo
     145
     146    write(*,*) "Downloading data for tsurf_slope done"
     147
     148#ifndef CPP_STD
     149    if (soil_pem) then
     150        write(*,*) "Downloading data for soiltemp_slope ..."
     151
     152        do islope = 1,nslope
     153            write(num,'(i2.2)') islope
     154            call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
     155        enddo
     156
     157        write(*,*) "Downloading data for soiltemp_slope done"
     158        write(*,*) "Downloading data for watersoil_density ..."
     159
     160        do islope = 1,nslope
     161            write(num,'(i2.2)') islope
     162            call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
     163        enddo
     164
     165        write(*,*) "Downloading data for  watersoil_density  done"
     166        write(*,*) "Downloading data for  watersurf_density  ..."
     167
     168        do islope = 1,nslope
     169            write(num,'(i2.2)') islope
     170            call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
     171        enddo
     172
     173        write(*,*) "Downloading data for  watersurf_density  done"
     174
     175    endif !soil_pem
     176#endif
     177else !nslope = 1 no slope, we copy all the values
     178    call get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
     179    call get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
     180    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
    135181   
    136  write(*,*) "Downloading data for tsurf_slope ..."
    137 
    138 DO islope=1,nslope
    139   write(num,fmt='(i2.2)') islope
    140   call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
    141 ENDDO
    142 
    143      write(*,*) "Downloading data for tsurf_slope done"
    144 
    145 #ifndef CPP_STD
    146 
    147      if(soil_pem) then
    148 
    149      write(*,*) "Downloading data for soiltemp_slope ..."
    150 
    151 DO islope=1,nslope
    152   write(num,fmt='(i2.2)') islope
    153   call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
    154 ENDDO
    155 
    156      write(*,*) "Downloading data for soiltemp_slope done"
    157 
    158      write(*,*) "Downloading data for watersoil_density ..."
    159 
    160 DO islope=1,nslope
    161   write(num,fmt='(i2.2)') islope
    162   call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
    163 ENDDO
    164 
    165      write(*,*) "Downloading data for  watersoil_density  done"
    166 
    167      write(*,*) "Downloading data for  watersurf_density  ..."
    168 
    169 DO islope=1,nslope
    170   write(num,fmt='(i2.2)') islope
    171   call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
    172 ENDDO
    173 
    174      write(*,*) "Downloading data for  watersurf_density  done"
    175 
    176   endif !soil_pem
    177 
    178 #endif
    179 
    180   else !nslope=1 no slope, we copy all the values
    181 
    182     CALL get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
    183     CALL get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
    184     call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
    185182#ifndef CPP_STD
    186183!    call get_var3("watercap", watercap_slope(:,:,1,:))
    187   watercap_slope(:,:,1,:)=0.
    188   watersurf_density_dyn(:,:,:,:)=0.
    189   watersoil_density_dyn(:,:,:,:,:)=0.
    190 #endif
    191 
    192     if(soil_pem) then
    193       call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
     184    watercap_slope(:,:,1,:) = 0.
     185    watersurf_density_dyn(:,:,:,:) = 0.
     186    watersoil_density_dyn(:,:,:,:,:) = 0.
     187#endif
     188
     189    if (soil_pem) call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
     190endif !nslope=1
     191
     192! Compute the minimum over the year for each point
     193write(*,*) "Computing the min of h2o_ice_slope"
     194min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap_slope,4)
     195write(*,*) "Computing the min of co2_ice_slope"
     196min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn,4)
     197
     198!Compute averages
     199write(*,*) "Computing average of tsurf"
     200tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
     201
     202#ifndef CPP_1D
     203    do islope = 1,nslope
     204        do t = 1,timelen
     205            call gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
     206        enddo
     207    enddo
     208#endif
     209
     210if (soil_pem) then
     211    write(*,*) "Computing average of tsoil"
     212    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
     213    write(*,*) "Computing average of watersurf_density"
     214    watersurf_density_ave(:,:) = sum(watersurf_density(:,:,:),3)/timelen
     215endif
     216
     217! By definition, a density is positive, we get rid of the negative values
     218where (min_co2_ice_dyn < 0.) min_co2_ice_dyn = 0.
     219where (min_h2o_ice_dyn < 0.) min_h2o_ice_dyn = 0.
     220
     221do i = 1,iim + 1
     222    do j = 1,jjm_input + 1
     223        do t = 1, timelen
     224            if (q_co2_dyn(i,j,t) < 0) then
     225                q_co2_dyn(i,j,t) = 1.e-10
     226            else if (q_co2_dyn(i,j,t) > 1) then
     227                q_co2_dyn(i,j,t) = 1.
     228            endif
     229            if (q_h2o_dyn(i,j,t) < 0) then
     230                q_h2o_dyn(i,j,t) = 1.e-30
     231            else if (q_h2o_dyn(i,j,t) > 1) then
     232                q_h2o_dyn(i,j,t) = 1.
     233            endif
     234            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
     235            vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
     236        enddo
     237    enddo
     238enddo
     239
     240#ifndef CPP_1D
     241    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
     242    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_GCM,ps_timeseries)
     243    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
     244    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
     245
     246    do islope = 1,nslope
     247        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
     248        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
     249        if (soil_pem) then
     250            do l = 1,nsoilmx
     251                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
     252                do t=1,timelen
     253                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
     254                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
     255                enddo
     256            enddo
     257        endif !soil_pem
     258        do t = 1,timelen
     259            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
     260            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
     261        enddo
     262    enddo
     263
     264    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_ave_dyn,tsurf_ave)
     265#else
     266    vmr_co2_gcm_phys(1,:) = vmr_co2_gcm(1,1,:)
     267    ps_timeseries(1,:) = ps_GCM(1,1,:)
     268    q_co2(1,:) = q_co2_dyn(1,1,:)
     269    q_h2o(1,:) = q_h2o_dyn(1,1,:)
     270    min_co2_ice(1,:) = min_co2_ice_dyn(1,1,:)
     271    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
     272    if (soil_pem) then
     273        tsoil_ave(1,:,:) = tsoil_ave_dyn(1,1,:,:)
     274        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
     275        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
    194276    endif !soil_pem
    195   endif !nslope=1
    196 
    197 ! Compute the minimum over the year for each point
    198   write(*,*) "Computing the min of h2o_ice_slope"
    199   min_h2o_ice_dyn(:,:,:)=minval(h2o_ice_s_dyn+watercap_slope,4)
    200   write(*,*) "Computing the min of co2_ice_slope"
    201   min_co2_ice_dyn(:,:,:)=minval(co2_ice_slope_dyn,4)
    202 
    203 !Compute averages
    204 
    205     write(*,*) "Computing average of tsurf"
    206     tsurf_ave_dyn(:,:,:)=SUM(tsurf_gcm_dyn(:,:,:,:),4)/timelen
    207 
    208 #ifndef CPP_1D
    209   DO islope = 1,nslope
    210     DO t=1,timelen
    211       CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
    212     ENDDO
    213   ENDDO
    214 #endif
    215 
    216   if(soil_pem) then
    217     write(*,*) "Computing average of tsoil"
    218     tsoil_ave_dyn(:,:,:,:)=SUM(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
    219     write(*,*) "Computing average of watersurf_density"
    220     watersurf_density_ave(:,:) = SUM(watersurf_density(:,:,:),3)/timelen
    221   endif
    222 
    223 ! By definition, a density is positive, we get rid of the negative values
    224   DO i=1,iim+1
    225     DO j = 1, jjm_input+1
    226        DO islope=1,nslope
    227           if (min_co2_ice_dyn(i,j,islope).LT.0) then
    228             min_co2_ice_dyn(i,j,islope)  = 0.
    229           endif
    230           if (min_h2o_ice_dyn(i,j,islope).LT.0) then
    231             min_h2o_ice_dyn(i,j,islope)  = 0.
    232           endif
    233        ENDDO
    234     ENDDO
    235   ENDDO
    236 
    237   DO i=1,iim+1
    238     DO j = 1, jjm_input+1
    239       DO t = 1, timelen
    240          if (q_co2_dyn(i,j,t).LT.0) then
    241               q_co2_dyn(i,j,t)=1E-10
    242          elseif (q_co2_dyn(i,j,t).GT.1) then
    243               q_co2_dyn(i,j,t)=1.
    244          endif
    245          if (q_h2o_dyn(i,j,t).LT.0) then
    246               q_h2o_dyn(i,j,t)=1E-30
    247          elseif (q_h2o_dyn(i,j,t).GT.1) then
    248               q_h2o_dyn(i,j,t)=1.
    249          endif
    250          mmean=1/(A*q_co2_dyn(i,j,t) +B)
    251          vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
    252       ENDDO
    253     ENDDO
    254   ENDDO
    255 
    256 #ifndef CPP_1D
    257 
    258      CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
    259      call gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,ps_GCM,ps_timeseries)
    260      CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,q_co2_dyn,q_co2)
    261      CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,q_h2o_dyn,q_h2o)
    262 
    263      DO islope = 1,nslope
    264        CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
    265        CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
    266        if(soil_pem) then
    267        DO l=1,nsoilmx
    268          CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
    269          DO t=1,timelen
    270            CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
    271            CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
    272          ENDDO
    273        ENDDO
    274        endif !soil_pem
    275        DO t=1,timelen
    276          CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
    277          CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
    278        ENDDO
    279      ENDDO
    280 
    281      CALL gr_dyn_fi(nslope,iim_input+1,jjm_input+1,ngrid,tsurf_ave_dyn,tsurf_ave)
    282 
    283 #else
    284 
    285   vmr_co2_gcm_phys(1,:)=vmr_co2_gcm(1,1,:)
    286   ps_timeseries(1,:)=ps_GCM(1,1,:)
    287   q_co2(1,:)=q_co2_dyn(1,1,:)
    288   q_h2o(1,:)=q_h2o_dyn(1,1,:)
    289   min_co2_ice(1,:)=min_co2_ice_dyn(1,1,:)
    290   min_h2o_ice(1,:)=min_h2o_ice_dyn(1,1,:)
    291   if(soil_pem) then
    292     tsoil_ave(1,:,:)=tsoil_ave_dyn(1,1,:,:)
    293     tsoil_gcm(1,:,:,:)=tsoil_gcm_dyn(1,1,:,:,:)
    294     watersoil_density(1,:,:,:)=watersoil_density_dyn(1,1,:,:,:)
    295   endif !soil_pem
    296   tsurf_GCM(1,:,:)=tsurf_GCM_dyn(1,1,:,:)
    297   co2_ice_slope(1,:,:)=co2_ice_slope_dyn(1,1,:,:)
    298   tsurf_ave(1,:)=tsurf_ave_dyn(1,1,:)
    299 
    300 #endif
    301 
    302   CONTAINS
     277    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
     278    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
     279    tsurf_ave(1,:) = tsurf_ave_dyn(1,1,:)
     280#endif
     281
     282END SUBROUTINE read_data_gcm
    303283
    304284SUBROUTINE check_dim(n1,n2,str1,str2)
    305   INTEGER,          INTENT(IN) :: n1, n2
    306   CHARACTER(LEN=*), INTENT(IN) :: str1, str2
    307   CHARACTER(LEN=256) :: s1, s2
    308   IF(n1/=n2) THEN
    309     s1='value of '//TRIM(str1)//' ='
    310     s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
    311     WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
    312     CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
    313   END IF
     285
     286implicit none
     287
     288integer,            intent(in) :: n1, n2
     289character(len = *), intent(in) :: str1, str2
     290
     291character(256) :: s1, s2
     292
     293if (n1 /= n2) then
     294    s1 = 'value of '//trim(str1)//' ='
     295    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
     296    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
     297    call abort_gcm(trim(modname),trim(msg),1)
     298endif
     299
    314300END SUBROUTINE check_dim
    315301
     302!=======================================================================
    316303
    317304SUBROUTINE get_var1(var,v)
    318   CHARACTER(LEN=*), INTENT(IN)  :: var
    319   REAL,             INTENT(OUT) :: v(:)
    320   CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
    321   CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
     305
     306implicit none
     307
     308character(len = *), intent(in)  :: var
     309real, dimension(:), intent(out) :: v
     310
     311call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
     312call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
     313
    322314END SUBROUTINE get_var1
    323315
     316!=======================================================================
    324317
    325318SUBROUTINE get_var3(var,v) ! on U grid
    326   CHARACTER(LEN=*), INTENT(IN)  :: var
    327   REAL,             INTENT(OUT) :: v(:,:,:)
    328   CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
    329   CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
     319
     320implicit none
     321
     322character(len = *),     intent(in)  :: var
     323real, dimension(:,:,:), intent(out) :: v
     324
     325call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
     326call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
    330327
    331328END SUBROUTINE get_var3
    332329
    333 SUBROUTINE get_var4(var,v)
    334   CHARACTER(LEN=*), INTENT(IN)  :: var
    335   REAL,             INTENT(OUT) :: v(:,:,:,:)
    336   CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
    337   CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
     330!=======================================================================
     331
     332SUBROUTINE get_var4(var,v)
     333
     334implicit none
     335
     336character(len = *),       intent(in)  :: var
     337real, dimension(:,:,:,:), intent(out) :: v
     338
     339call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
     340call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
     341
    338342END SUBROUTINE get_var4
    339343
    340 SUBROUTINE err(ierr,typ,nam)
    341   INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
    342   CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
    343   CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
    344   IF(ierr==NF90_NoERR) RETURN
    345   SELECT CASE(typ)
    346     CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
    347     CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
    348     CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
    349     CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
    350   END SELECT
    351   CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
    352 END SUBROUTINE err
    353 
    354 END SUBROUTINE read_data_gcm
     344!=======================================================================
     345
     346SUBROUTINE error_msg(ierr,typ,nam)
     347
     348implicit none
     349
     350integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
     351character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
     352character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
     353
     354if (ierr == nf90_noerr) return
     355select case(typ)
     356    case('inq');   msg="Field <"//trim(nam)//"> is missing"
     357    case('get');   msg="Reading failed for <"//trim(nam)//">"
     358    case('open');  msg="File opening failed for <"//trim(nam)//">"
     359    case('close'); msg="File closing failed for <"//trim(nam)//">"
     360    case default
     361        write(*,*) 'There is no message for this error.'
     362        error stop
     363end select
     364call abort_gcm(trim(modname),trim(msg),ierr)
     365
     366END SUBROUTINE error_msg
     367
     368END MODULE read_data_GCM_mod

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

    r3069 r3076  
    66! Author: RV, JBC
    77!=======================================================================
    8 IMPLICIT NONE
    98
    10 CONTAINS
     9implicit none
     10
     11!=======================================================================
     12contains
     13!=======================================================================
    1114
    1215SUBROUTINE recomp_orb_param(i_myear,year_iter)
    1316
    14   use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
    15 #ifndef CPP_STD     
     17use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
     18use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, end_lask_param_mod, last_ilask
     19#ifndef CPP_STD
    1620    use comconst_mod, only: pi
    1721    use planete_h,    only: e_elips, obliquit, lsperi, periheli, aphelie, p_elips, peri_day, year_day
    … …  
    2125#endif
    2226
    23   use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, end_lask_param_mod, last_ilask
    24 
    25   IMPLICIT NONE
     27implicit none
    2628
    2729!--------------------------------------------------------
    2830! Input Variables
    2931!--------------------------------------------------------
    30   integer, intent(in) :: i_myear   ! Number of simulated Martian years
    31   integer, intent(in) :: year_iter ! Number of iterations of the PEM
     32integer, intent(in) :: i_myear   ! Number of simulated Martian years
     33integer, intent(in) :: year_iter ! Number of iterations of the PEM
    3234
    3335!--------------------------------------------------------
    … …  
    3638
    3739!--------------------------------------------------------
    38 ! Local variables 
     40! Local variables
    3941!--------------------------------------------------------
    40   integer :: Year           ! Year of the simulation
    41   real    :: Lsp            ! Ls perihelion
    42   integer :: ilask          ! Loop variable
    43   real    :: halfaxe        ! Million of km
    44   real    :: unitastr
    45   real    :: xa, xb, ya, yb ! Variables for interpolation
    46   logical :: found_year     ! Flag variable
     42integer :: Year           ! Year of the simulation
     43real    :: Lsp            ! Ls perihelion
     44integer :: ilask          ! Loop variable
     45real    :: halfaxe        ! Million of km
     46real    :: unitastr
     47real    :: xa, xb, ya, yb ! Variables for interpolation
     48logical :: found_year     ! Flag variable
    4749
    4850! **********************************************************************
    4951! 0. Initializations
    5052! **********************************************************************
    51 #ifdef CPP_STD
    52     real aphelie
    53     real periheli
     53#ifdef CPP_STD
     54    real :: aphelie, periheli
    5455
    55     aphelie=apoastr
    56     periheli=periastr
     56    aphelie = apoastr
     57    periheli = periastr
    5758#endif
    5859
    … …  
    111112p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
    112113
    113 #ifndef CPP_STD 
     114#ifndef CPP_STD
    114115    call call_dayperi(Lsp,e_elips,peri_day,year_day)
    115116#endif
    … …  
    122123END SUBROUTINE recomp_orb_param
    123124
    124 !********************************************************************************   
    125      
     125!********************************************************************************
     126
    126127END MODULE recomp_orb_param_mod
    127128

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

    r3075 r3076  
    1 !
    2 ! $Id $
    3 !
    4 SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope, emissivity_slope, &
     1MODULE recomp_tend_co2_slope_mod
     2
     3implicit none
     4
     5!=======================================================================
     6contains
     7!=======================================================================
     8
     9SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope,emissivity_slope, &
    510                                 vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new)
    611
    7       use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB,Lco2, sols_per_my, sec_per_sol
     12use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB,Lco2, sols_per_my, sec_per_sol
    813
    9       IMPLICIT NONE
     14implicit none
    1015
    1116!=======================================================================
    … …  
    1722!   arguments:
    1823!   ----------
    19 
    2024!   INPUT
    21   INTEGER, intent(in) :: timelen,ngrid,nslope
    22   REAL, INTENT(in) ::  vmr_co2_gcm(ngrid,timelen)                ! physical point field : Volume mixing ratio of co2 in the first layer
    23   REAL, INTENT(in) ::  vmr_co2_pem(ngrid,timelen)                ! physical point field : Volume mixing ratio of co2 in the first layer
    24   REAL, intent(in) :: ps_GCM_2(ngrid,timelen)                    ! physical point field : Surface pressure in the GCM
    25   REAL, intent(in) :: global_ave_press_GCM                       ! global averaged pressure at previous timestep
    26   REAL, intent(in) :: global_ave_press_new                       ! global averaged pressure at current timestep
    27   REAL, intent(in) ::  tendencies_co2_ice_phys_ini(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    28   REAL, intent(in) :: co2ice_slope(ngrid,nslope)                 ! CO2 ice per mesh and sub-grid slope(kg/m^2)
    29   REAL, intent(in) :: emissivity_slope(ngrid,nslope)             ! Emissivity per mesh and sub-grid slope(1)
    30 
     25integer,                        intent(in) :: timelen, ngrid, nslope
     26real, dimension(ngrid,timelen), intent(in) :: vmr_co2_gcm                 ! physical point field : Volume mixing ratio of co2 in the first layer
     27real, dimension(ngrid,timelen), intent(in) :: vmr_co2_pem                 ! physical point field : Volume mixing ratio of co2 in the first layer
     28real, dimension(ngrid,timelen), intent(in) :: ps_GCM_2                    ! physical point field : Surface pressure in the GCM
     29real,                           intent(in) :: global_ave_press_GCM        ! global averaged pressure at previous timestep
     30real,                           intent(in) :: global_ave_press_new        ! global averaged pressure at current timestep
     31real, dimension(ngrid,nslope),  intent(in) :: tendencies_co2_ice_phys_ini ! physical point field : Evolution of perenial ice over one year
     32real, dimension(ngrid,nslope),  intent(in) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope(kg/m^2)
     33real, dimension(ngrid,nslope),  intent(in) :: emissivity_slope            ! Emissivity per mesh and sub-grid slope(1)
    3134!   OUTPUT
    32   REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
     35real, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    3336
    3437!   local:
    3538!   ----
    36 
    37   INTEGER :: i,t,islope
    38   REAL ::  coef, ave
     39integer :: i, t, islope
     40real    :: coef, ave
    3941
    4042! Evolution of the water ice for each physical point
    41   do i=1,ngrid
    42     do islope=1,nslope
    43       coef=sols_per_my*sec_per_sol*emissivity_slope(i,islope)*sigmaB/Lco2
    44       ave=0.
    45       if(co2ice_slope(i,islope).gt.1e-4 .and. abs(tendencies_co2_ice_phys(i,islope)).gt.1e-5) then
    46         do t=1,timelen
    47            ave=ave+(beta_clap_co2/(alpha_clap_co2-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4 &
    48               -(beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*(global_ave_press_new/global_ave_press_GCM)/100.)))**4
    49          enddo
    50       endif
    51       if(ave.lt.1e-4) ave = 0.
    52       tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)-coef*ave/timelen
     43do i = 1,ngrid
     44    do islope = 1,nslope
     45        coef = sols_per_my*sec_per_sol*emissivity_slope(i,islope)*sigmaB/Lco2
     46        ave = 0.
     47        if (co2ice_slope(i,islope) > 1.e-4 .and. abs(tendencies_co2_ice_phys(i,islope)) > 1.e-5) then
     48            do t=1,timelen
     49                ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4 &
     50                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*(global_ave_press_new/global_ave_press_GCM)/100.)))**4
     51            enddo
     52        endif
     53        if (ave < 1e-4) ave = 0.
     54        tendencies_co2_ice_phys(i,islope) = tendencies_co2_ice_phys_ini(i,islope) - coef*ave/timelen
    5355    enddo
    54   enddo
     56enddo
    5557
    5658END SUBROUTINE recomp_tend_co2_slope
     59
     60END MODULE recomp_tend_co2_slope_mod

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

    r3050 r3076  
    1 program reshape_XIOS_output
     1PROGRAM reshape_XIOS_output
    22
    33!=======================================================================
    … …  
    1010! Authors: RV & LL
    1111!=======================================================================
    12     use netcdf
    13     implicit none
    14     integer :: status, ncid, ncid1, ncid2
    15     integer :: nDims, nVars, nGlobalAtts, unlimDimID
    16     integer i,j
    17 
    18     integer :: include_parents
    19 
    20     integer, dimension(:),allocatable :: dimids
    21     integer, dimension(:),allocatable :: varids
    22 
    23     integer, dimension(:),allocatable :: dimids_2
    24     integer, dimension(:),allocatable :: varids_2
    25 
    26     integer, dimension(:),allocatable :: dimid_var
    27 
    28     real, dimension(:), allocatable :: tempvalues_1d
    29     real, dimension(:), allocatable :: values_1d
    30 
    31     real, dimension(:,:), allocatable :: tempvalues_2d
    32     real, dimension(:,:), allocatable :: values_2d
    33 
    34     real, dimension(:,:,:), allocatable :: tempvalues_3d
    35     real, dimension(:,:,:), allocatable :: values_3d
    36 
    37     real, dimension(:,:,:,:), allocatable :: tempvalues_4d
    38     real, dimension(:,:,:,:), allocatable :: values_4d
    39 
    40   character*1 str2
    41   character*30 :: name_
    42   character*30 :: namevar
    43   integer  :: xtype_var
    44   integer :: len_
    45   integer :: len_1,len_2
    46   integer :: len_lat, len_lon, len_time, len_soil
    47   integer :: dimid_lon, dimid_lat, dimid_time, dimid_soil
    48   integer :: dimid_2
    49   integer :: numdims
    50   integer :: numatts
    51   integer :: numyear
    52 
    53 DO numyear=1, 2
    54 write(*,*) 'numyear',numyear
    55 write(str2(1:1),'(i1.1)') numyear
    56 !nf90_open                 ! open existing netCDF dataset
    57 !integer :: ncid, status
    58 !...
    59 status = nf90_open(path = "data2reshape"//str2//".nc", mode = nf90_nowrite, ncid = ncid1)
    60 if(status /= nf90_noerr) call handle_err(status)
    61 
    62 status = nf90_create(path = "datareshaped"//str2//".nc", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2)
    63 if(status /= nf90_noerr) call handle_err(status)
    64 
    65 status = nf90_inquire(ncid1, ndims, nvars, nglobalatts, unlimdimid)
    66 if(status /= nf90_noerr) call handle_err(status)
    67 
    68 allocate(dimids(ndims))
    69 allocate(varids(nvars))
    70 
    71 allocate(dimids_2(ndims))
    72 allocate(varids_2(nvars))
    73 
    74 status = nf90_inq_dimids(ncid1, ndims, dimids, include_parents)
    75 if(status /= nf90_noerr) call handle_err(status)
    76 status = nf90_inq_varids(ncid1, nvars, varids)
    77 if(status /= nf90_noerr) call handle_err(status)
    78 
    79 do i=1,ndims
    80   status = nf90_inquire_dimension(ncid1, dimids(i), name_, len_)
    81   if(status /= nf90_noerr) call handle_err(status)
    82   if(name_.eq."lon" .or. name_.eq."longitude")  then
    83      dimid_lon=dimids(i)
    84      len_lon=len_
    85      len_=len_+1
    86   elseif(name_.eq."lat".or. name_.eq."latitude") then
    87      dimid_lat=dimids(i)
    88      len_lat=len_
    89   elseif(name_.eq."time_counter".or. name_.eq. "Time") then
    90      dimid_time=dimids(i)
    91      len_time=len_
    92   elseif(name_.eq."soil_layers".or. name_.eq. "subsurface_layers") then
    93      dimid_soil=dimids(i)
    94      len_soil=len_
    95   endif
    96   status = nf90_def_dim(ncid2, name_, len_, dimid_2)
    97   if(status /= nf90_noerr) call handle_err(status)
    98   dimids_2(i)=dimid_2
     12use netcdf
     13
     14implicit none
     15
     16integer                               :: state, ncid, ncid1, ncid2, nDims, nVars, nGlobalAtts, unlimDimID
     17integer                               :: i, j, include_parents
     18integer, dimension(:),    allocatable :: dimids, varids, dimids_2, varids_2, dimid_var
     19real, dimension(:),       allocatable :: tempvalues_1d, values_1d
     20real, dimension(:,:),     allocatable :: tempvalues_2d, values_2d
     21real, dimension(:,:,:),   allocatable :: tempvalues_3d, values_3d
     22real, dimension(:,:,:,:), allocatable :: tempvalues_4d, values_4d
     23character(1)                          :: str2
     24character(30)                         :: name_, namevar
     25integer                               :: xtype_var, len_, len_1, len_2, len_lat, len_lon, len_time, len_soil
     26integer                               :: dimid_lon, dimid_lat, dimid_time, dimid_soil, dimid_2, numdims, numatts, numyear
     27
     28do numyear = 1,2
     29    write(*,*) 'numyear',numyear
     30    write(str2(1:1),'(i1.1)') numyear
     31    !nf90_open                 ! open existing netCDF dataset
     32    !integer :: ncid, state
     33    !...
     34    state = nf90_open(path = "data2reshape"//str2//".nc", mode = nf90_nowrite, ncid = ncid1)
     35    if (state /= nf90_noerr) call handle_err(state)
     36
     37    state = nf90_create(path = "datareshaped"//str2//".nc", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2)
     38    if (state /= nf90_noerr) call handle_err(state)
     39
     40    state = nf90_inquire(ncid1, ndims, nvars, nglobalatts, unlimdimid)
     41    if (state /= nf90_noerr) call handle_err(state)
     42
     43    allocate(dimids(ndims))
     44    allocate(varids(nvars))
     45
     46    allocate(dimids_2(ndims))
     47    allocate(varids_2(nvars))
     48
     49    state = nf90_inq_dimids(ncid1, ndims, dimids, include_parents)
     50    if (state /= nf90_noerr) call handle_err(state)
     51    state = nf90_inq_varids(ncid1, nvars, varids)
     52    if (state /= nf90_noerr) call handle_err(state)
     53
     54    do i = 1,ndims
     55        state = nf90_inquire_dimension(ncid1, dimids(i), name_, len_)
     56        if (state /= nf90_noerr) call handle_err(state)
     57        if (name_ == "lon" .or. name_ == "longitude") then
     58            dimid_lon = dimids(i)
     59            len_lon = len_
     60            len_ = len_ + 1
     61        elseif (name_ == "lat".or. name_ == "latitude") then
     62            dimid_lat=dimids(i)
     63            len_lat=len_
     64        elseif (name_ == "time_counter".or. name_ ==  "Time") then
     65            dimid_time=dimids(i)
     66            len_time=len_
     67        elseif (name_ == "soil_layers".or. name_ ==  "subsurface_layers") then
     68            dimid_soil=dimids(i)
     69            len_soil = len_
     70        endif
     71        state = nf90_def_dim(ncid2, name_,len_,dimid_2)
     72        if (state /= nf90_noerr) call handle_err(state)
     73        dimids_2(i) = dimid_2
     74    enddo
     75
     76    do i = 1,nvars
     77        state = nf90_inquire_variable(ncid1, varids(i),name = namevar,xtype = xtype_var,ndims = numdims,natts = numatts)
     78        write(*,*) "namevar00= ", namevar
     79        if (state /= nf90_noerr) call handle_err(state)
     80        allocate(dimid_var(numdims))
     81        state = nf90_inquire_variable(ncid1,varids(i),name = namevar,xtype = xtype_var,ndims = numdims,dimids = dimid_var,natts = numatts)
     82        if (state /= nf90_noerr) call handle_err(state)
     83        if (numdims == 1) then
     84            if (namevar == "lon") then
     85                allocate(tempvalues_1d(len_lon))
     86                allocate(values_1d(len_lon + 1))
     87                state = nf90_get_var(ncid1,varids(i),tempvalues_1d)
     88                if (state /= nf90_noerr) call handle_err(state)
     89                state = nf90_def_var(ncid2,namevar,xtype_var, dimid_var, varids_2(i))
     90                if (state /= nf90_noerr) call handle_err(state)
     91                values_1d(1:len_lon) = tempvalues_1d(:)
     92                values_1d(len_lon + 1) = values_1d(1)
     93                state = nf90_enddef(ncid2)
     94                if (state /= nf90_noerr) call handle_err(state)
     95                state = nf90_put_var(ncid2, varids_2(i), values_1d)
     96                if (state /= nf90_noerr) call handle_err(state)
     97                state = nf90_redef(ncid2)
     98                if (state /= nf90_noerr) call handle_err(state)
     99                deallocate(tempvalues_1d)
     100                deallocate(values_1d)
     101            else
     102                state = nf90_inquire_dimension(ncid1,dimid_var(1),name_,len_)
     103                if (state /= nf90_noerr) call handle_err(state)
     104                allocate(tempvalues_1d(len_))
     105                state = nf90_get_var(ncid1,varids(i),tempvalues_1d)
     106                if (state /= nf90_noerr) call handle_err(state)
     107                state = nf90_def_var(ncid2,namevar,xtype_var,dimid_var,varids_2(i))
     108                if (state /= nf90_noerr) call handle_err(state)
     109                state = nf90_enddef(ncid2)
     110                if (state /= nf90_noerr) call handle_err(state)
     111                state = nf90_put_var(ncid2, varids_2(i), tempvalues_1d)
     112                if (state /= nf90_noerr) call handle_err(state)
     113                state = nf90_redef(ncid2)
     114                if (state /= nf90_noerr) call handle_err(state)
     115                deallocate(tempvalues_1d)
     116            endif
     117        else if (numdims == 2) then
     118            if (namevar == "area") then
     119                allocate(tempvalues_2d(len_lon,len_lat))
     120                allocate(values_2d(len_lon + 1,len_lat))
     121                state = nf90_get_var(ncid1,varids(i),tempvalues_2d)
     122                if (state /= nf90_noerr) call handle_err(state)
     123                state = nf90_def_var(ncid2,namevar,xtype_var,dimid_var,varids_2(i))
     124                if (state /= nf90_noerr) call handle_err(state)
     125                values_2d(1:len_lon,:) = tempvalues_2d(:,:)
     126                values_2d(len_lon+1,:) = values_2d(1,:)
     127                state = nf90_enddef(ncid2)
     128                if (state /= nf90_noerr) call handle_err(state)
     129                state = nf90_put_var(ncid2,varids_2(i),values_2d)
     130                if (state /= nf90_noerr) call handle_err(state)
     131                state = nf90_redef(ncid2)
     132                if (state /= nf90_noerr) call handle_err(state)
     133                deallocate(tempvalues_2d)
     134                deallocate(values_2d)
     135            else
     136                state = nf90_inquire_dimension(ncid1,dimid_var(1),name_,len_1)
     137                if (state /= nf90_noerr) call handle_err(state)
     138                state = nf90_inquire_dimension(ncid1,dimid_var(2),name_,len_2)
     139                if (state /= nf90_noerr) call handle_err(state)
     140                allocate(tempvalues_2d(len_1,len_2))
     141                state = nf90_get_var(ncid1, varids(i), tempvalues_2d)
     142                if (state /= nf90_noerr) call handle_err(state)
     143                state = nf90_def_var(ncid2,namevar,xtype_var,dimid_var,varids_2(i))
     144                if (state /= nf90_noerr) call handle_err(state)
     145                state = nf90_enddef(ncid2)
     146                if (state /= nf90_noerr) call handle_err(state)
     147                state = nf90_put_var(ncid2, varids_2(i), tempvalues_2d)
     148                if (state /= nf90_noerr) call handle_err(state)
     149                state = nf90_redef(ncid2)
     150                if (state /= nf90_noerr) call handle_err(state)
     151                deallocate(tempvalues_2d)
     152            endif
     153        elseif (numdims == 3) then
     154            allocate(tempvalues_3d(len_lon,len_lat,len_time))
     155            allocate(values_3d(len_lon + 1,len_lat,len_time))
     156            state = nf90_get_var(ncid1,varids(i),tempvalues_3d)
     157            if (state /= nf90_noerr) call handle_err(state)
     158            state = nf90_def_var(ncid2,namevar,xtype_var,dimid_var,varids_2(i))
     159            if (state /= nf90_noerr) call handle_err(state)
     160            values_3d(1:len_lon,:,:) = tempvalues_3d(:,:,:)
     161            values_3d(len_lon+1,:,:) = values_3d(1,:,:)
     162            state = nf90_enddef(ncid2)
     163            if (state /= nf90_noerr) call handle_err(state)
     164            state = nf90_put_var(ncid2, varids_2(i), values_3d)
     165            if (state /= nf90_noerr) call handle_err(state)
     166            state = nf90_redef(ncid2)
     167            if (state /= nf90_noerr) call handle_err(state)
     168            deallocate(tempvalues_3d)
     169            deallocate(values_3d)
     170        else if (numdims == 4) then
     171            allocate(tempvalues_4d(len_lon,len_lat,len_soil,len_time))
     172            allocate(values_4d(len_lon+1,len_lat,len_soil,len_time))
     173            state = nf90_get_var(ncid1, varids(i), tempvalues_4d)
     174            if (state /= nf90_noerr) call handle_err(state)
     175            state = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
     176            if (state /= nf90_noerr) call handle_err(state)
     177            state = nf90_enddef(ncid2)
     178            values_4d(1:len_lon,:,:,:) = tempvalues_4d(:,:,:,:)
     179            values_4d(len_lon+1,:,:,:) = values_4d(1,:,:,:)
     180            if (state /= nf90_noerr) call handle_err(state)
     181            state = nf90_put_var(ncid2, varids_2(i), values_4d)
     182            if (state /= nf90_noerr) call handle_err(state)
     183            state = nf90_redef(ncid2)
     184            if (state /= nf90_noerr) call handle_err(state)
     185            deallocate(tempvalues_4d)
     186            deallocate(values_4d)
     187        endif
     188        deallocate(dimid_var)
     189    enddo
     190
     191    state = nf90_enddef(ncid2)
     192    if (state /= nf90_noerr) call handle_err(state)
     193    state = nf90_close(ncid1)
     194    if (state /= nf90_noerr) call handle_err(state)
     195    state = nf90_close(ncid2)
     196    if (state /= nf90_noerr) call handle_err(state)
     197
     198    deallocate(dimids,varids,dimids_2,varids_2)
    99199enddo
    100200
    101 do i=1,nvars
    102   status = nf90_inquire_variable(ncid1, varids(i), name=namevar, xtype=xtype_var, ndims = numdims,natts = numatts)
    103   write(*,*) "namevar00= ", namevar
    104   if(status /= nf90_noerr) call handle_err(status)
    105   allocate(dimid_var(numdims))
    106   status = nf90_inquire_variable(ncid1, varids(i), name=namevar, xtype=xtype_var, ndims = numdims, dimids=dimid_var, natts = numatts)
    107   if(status /= nf90_noerr) call handle_err(status)
    108   if(numdims.eq.1) then
    109     if(namevar.eq."lon") then
    110       allocate(tempvalues_1d(len_lon))
    111       allocate(values_1d(len_lon+1))
    112       status = nf90_get_var(ncid1, varids(i), tempvalues_1d)
    113       if(status /= nf90_noerr) call handle_err(status)
    114       status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
    115       if(status /= nf90_noerr) call handle_err(status)
    116       values_1d(1:len_lon)=tempvalues_1d(:)
    117       values_1d(len_lon+1)=values_1d(1)
    118       status = nf90_enddef(ncid2)
    119       if(status /= nf90_noerr) call handle_err(status)
    120       status = nf90_put_var(ncid2, varids_2(i), values_1d)
    121       if(status /= nf90_noerr) call handle_err(status)
    122       status = nf90_redef(ncid2)
    123       if(status /= nf90_noerr) call handle_err(status)
    124       deallocate(tempvalues_1d)
    125       deallocate(values_1d) 
    126     else
    127       status = nf90_inquire_dimension(ncid1, dimid_var(1), name_, len_)
    128       if(status /= nf90_noerr) call handle_err(status)
    129       allocate(tempvalues_1d(len_))
    130       status = nf90_get_var(ncid1, varids(i), tempvalues_1d)
    131       if(status /= nf90_noerr) call handle_err(status)
    132       status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
    133       if(status /= nf90_noerr) call handle_err(status)
    134       status = nf90_enddef(ncid2)
    135       if(status /= nf90_noerr) call handle_err(status)
    136       status = nf90_put_var(ncid2, varids_2(i), tempvalues_1d) 
    137       if(status /= nf90_noerr) call handle_err(status)
    138       status = nf90_redef(ncid2)
    139       if(status /= nf90_noerr) call handle_err(status)
    140       deallocate(tempvalues_1d)   
    141     endif
    142   elseif(numdims.eq.2) then
    143     if(namevar.eq."area") then
    144       allocate(tempvalues_2d(len_lon,len_lat))
    145       allocate(values_2d(len_lon+1,len_lat))     
    146       status = nf90_get_var(ncid1, varids(i), tempvalues_2d)
    147       if(status /= nf90_noerr) call handle_err(status)
    148       status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
    149       if(status /= nf90_noerr) call handle_err(status)
    150       values_2d(1:len_lon,:)=tempvalues_2d(:,:)
    151       values_2d(len_lon+1,:)=values_2d(1,:)
    152       status = nf90_enddef(ncid2)
    153       if(status /= nf90_noerr) call handle_err(status)
    154       status = nf90_put_var(ncid2, varids_2(i), values_2d)   
    155       if(status /= nf90_noerr) call handle_err(status)
    156       status = nf90_redef(ncid2)
    157       if(status /= nf90_noerr) call handle_err(status)
    158       deallocate(tempvalues_2d)
    159       deallocate(values_2d)
    160     else
    161       status = nf90_inquire_dimension(ncid1, dimid_var(1), name_, len_1)
    162       if(status /= nf90_noerr) call handle_err(status)
    163       status = nf90_inquire_dimension(ncid1, dimid_var(2), name_, len_2)
    164       if(status /= nf90_noerr) call handle_err(status)
    165       allocate(tempvalues_2d(len_1,len_2))
    166       status = nf90_get_var(ncid1, varids(i), tempvalues_2d)
    167       if(status /= nf90_noerr) call handle_err(status)
    168       status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
    169       if(status /= nf90_noerr) call handle_err(status)
    170       status = nf90_enddef(ncid2)
    171       if(status /= nf90_noerr) call handle_err(status)
    172       status = nf90_put_var(ncid2, varids_2(i), tempvalues_2d)
    173       if(status /= nf90_noerr) call handle_err(status)
    174       status = nf90_redef(ncid2)
    175       if(status /= nf90_noerr) call handle_err(status)
    176       deallocate(tempvalues_2d)
    177     endif
    178   elseif(numdims.eq.3) then
    179       allocate(tempvalues_3d(len_lon,len_lat,len_time))
    180       allocate(values_3d(len_lon+1,len_lat,len_time))
    181       status = nf90_get_var(ncid1, varids(i), tempvalues_3d)
    182       if(status /= nf90_noerr) call handle_err(status)
    183       status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
    184       if(status /= nf90_noerr) call handle_err(status)
    185       values_3d(1:len_lon,:,:)=tempvalues_3d(:,:,:)
    186       values_3d(len_lon+1,:,:)=values_3d(1,:,:)
    187       status = nf90_enddef(ncid2)
    188       if(status /= nf90_noerr) call handle_err(status)
    189       status = nf90_put_var(ncid2, varids_2(i), values_3d)
    190       if(status /= nf90_noerr) call handle_err(status)
    191       status = nf90_redef(ncid2)
    192       if(status /= nf90_noerr) call handle_err(status)
    193       deallocate(tempvalues_3d)
    194       deallocate(values_3d)
    195   elseif(numdims.eq.4) then
    196       allocate(tempvalues_4d(len_lon,len_lat,len_soil,len_time))
    197       allocate(values_4d(len_lon+1,len_lat,len_soil,len_time))
    198       status = nf90_get_var(ncid1, varids(i), tempvalues_4d)
    199       if(status /= nf90_noerr) call handle_err(status)
    200       status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
    201       if(status /= nf90_noerr) call handle_err(status)
    202       status = nf90_enddef(ncid2)
    203       values_4d(1:len_lon,:,:,:)=tempvalues_4d(:,:,:,:)
    204       values_4d(len_lon+1,:,:,:)=values_4d(1,:,:,:)
    205       if(status /= nf90_noerr) call handle_err(status)
    206       status = nf90_put_var(ncid2, varids_2(i), values_4d)
    207       if(status /= nf90_noerr) call handle_err(status)
    208       status = nf90_redef(ncid2)
    209       if(status /= nf90_noerr) call handle_err(status)
    210       deallocate(tempvalues_4d)
    211       deallocate(values_4d)
    212   endif
    213 
    214   deallocate(dimid_var)
    215 enddo
    216 
    217 status = nf90_enddef(ncid2)
    218 if(status /= nf90_noerr) call handle_err(status)
    219 status = nf90_close(ncid1)
    220 if(status /= nf90_noerr) call handle_err(status)
    221 status = nf90_close(ncid2)
    222 if(status /= nf90_noerr) call handle_err(status)
    223 
    224 
    225 deallocate(dimids)
    226 deallocate(varids)
    227 deallocate(dimids_2)
    228 deallocate(varids_2)
    229 
    230 enddo
    231 
    232 end program reshape_XIOS_output
    233 
     201END PROGRAM reshape_XIOS_output
     202

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

    r3075 r3076  
    1       SUBROUTINE soil_TIfeedback_PEM(ngrid,nsoil,icecover,   newtherm_i)
     1MODULE soil_TIfeedback_PEM_mod
    22
    3       use comsoil_h_PEM, only: layer_PEM, inertiedat_PEM
     3implicit none
    44
    5       IMPLICIT NONE
     5!=======================================================================
     6contains
     7!=======================================================================
     8
     9SUBROUTINE soil_TIfeedback_PEM(ngrid,nsoil,icecover,newtherm_i)
     10
     11use comsoil_h_PEM, only: layer_PEM, inertiedat_PEM
     12
     13implicit none
    614
    715!=======================================================================
    … …  
    1321!   - One layer of surface water ice (the thickness is given
    1422!     by the variable icecover (in kg of ice per m2) and the thermal
    15 !     inertia is prescribed by inert_h2o_ice (see surfdat_h)); 
     23!     inertia is prescribed by inert_h2o_ice (see surfdat_h));
    1624!   - A transitional layer of mixed thermal inertia;
    1725!   - A last layer of regolith below the ice cover whose thermal inertia
    … …  
    2432!=======================================================================
    2533
     34!Inputs
     35!------
     36integer,                intent(in) :: ngrid    ! Number of horizontal grid points
     37integer,                intent(in) :: nsoil    ! Number of soil layers
     38real, dimension(ngrid), intent(in) :: icecover ! tracer on the surface (kg.m-2)
     39!Outputs
     40!-------
     41real,intent(inout) :: newtherm_i(ngrid,nsoil)    ! New soil thermal inertia
    2642!Local variables
    2743!---------------
     44integer :: ig                       ! Grid point (ngrid)
     45integer :: ik                       ! Grid point (nsoil)
     46integer :: iref                     ! Ice/Regolith boundary index
     47real    :: icedepth                 ! Ice cover thickness (m)
     48real    :: inert_h2o_ice = 800.     ! surface water ice thermal inertia [SI]
     49real    :: rho_ice = 920.           ! density of water ice [kg/m^3]
     50real    :: prev_thermi(ngrid,nsoil) ! previous thermal inertia [SI]
    2851
    29       INTEGER :: ig                     ! Grid point (ngrid)
    30       INTEGER :: ik                     ! Grid point (nsoil)
    31       INTEGER :: iref                   ! Ice/Regolith boundary index
    32       INTEGER, INTENT(IN) :: ngrid                  ! Number of horizontal grid points
    33       INTEGER, INTENT(IN) :: nsoil                  ! Number of soil layers
    34       REAL :: icedepth                  ! Ice cover thickness (m)
    35       REAL :: inert_h2o_ice = 800.      ! surface water ice thermal inertia [SI]
    36       REAL :: rho_ice = 920.            ! density of water ice [kg/m^3]
    37       REAL :: prev_thermi(ngrid,nsoil)  ! previous thermal inertia [SI]
    38 !Inputs
    39 !------
    40 
    41       REAL ,INTENT(IN):: icecover(ngrid)         ! tracer on the surface (kg.m-2)
    42 
    43 !Outputs
    44 !-------
    45 
    46       REAL,INTENT(INOUT) :: newtherm_i(ngrid,nsoil)    ! New soil thermal inertia
    47 
    48       prev_thermi(:,:) = newtherm_i(:,:)
     52prev_thermi(:,:) = newtherm_i(:,:)
    4953
    5054!Creating the new soil thermal inertia table
    5155!-------------------------------------------
    52       DO ig=1,ngrid
    53 !      Calculating the ice cover thickness
    54        
    55         icedepth=icecover(ig)/rho_ice
    56        
    57 !      If the ice cover is too thick or watercaptag=true,
    58 !        the entire column is changed :
    59         IF (icedepth.ge.layer_PEM(nsoil)) THEN
    60           DO ik=1,nsoil
    61                newtherm_i(ig,ik)=max(inert_h2o_ice,prev_thermi(ig,ik))
    62           ENDDO
    63 !      We neglect the effect of a very thin ice cover :
    64         ELSE IF (icedepth.lt.layer_PEM(1)) THEN
    65           DO ik=1,nsoil
    66                newtherm_i(ig,ik)=inertiedat_PEM(ig,ik)
    67           ENDDO
    68         ELSE
    69 !        Ice/regolith boundary index :
    70           iref=1
    71 !        Otherwise, we find the ice/regolith boundary:
    72           DO ik=1,nsoil-1
    73               IF ((icedepth.ge.layer_PEM(ik)).and. (icedepth.lt.layer_PEM(ik+1))) THEN
    74                   iref=ik+1
    75                   EXIT
    76               ENDIF
    77           ENDDO
    78 !        And we change the thermal inertia:
    79           DO ik=1,iref-1
    80             newtherm_i(ig,ik)=max(inert_h2o_ice,prev_thermi(ig,ik))
    81           ENDDO
    82 !        Transition (based on the equations of thermal conduction):
    83           newtherm_i(ig,iref)=sqrt( (layer_PEM(iref)-layer_PEM(iref-1)) / &
    84            ( ((icedepth-layer_PEM(iref-1))/newtherm_i(ig,iref-1)**2) + &
    85              ((layer_PEM(iref)-icedepth)/inertiedat_PEM(ig,ik)**2) ) )
    86 !        Underlying regolith:
    87           DO ik=iref+1,nsoil
    88               newtherm_i(ig,ik)=inertiedat_PEM(ig,ik)
    89           ENDDO
    90         ENDIF ! icedepth
    91       ENDDO ! ig
     56do ig = 1,ngrid
     57    ! Calculating the ice cover thickness
     58    icedepth = icecover(ig)/rho_ice
    9259
    93 !=======================================================================
    94       RETURN
    95       END
     60    ! If the ice cover is too thick or watercaptag = true, the entire column is changed:
     61    if (icedepth >= layer_PEM(nsoil)) then
     62        do ik = 1,nsoil
     63            newtherm_i(ig,ik) = max(inert_h2o_ice,prev_thermi(ig,ik))
     64        enddo
     65    ! We neglect the effect of a very thin ice cover:
     66    else if (icedepth < layer_PEM(1)) then
     67        do ik = 1,nsoil
     68            newtherm_i(ig,ik) = inertiedat_PEM(ig,ik)
     69        enddo
     70    else
     71        ! Ice/regolith boundary index:
     72        iref = 1
     73        ! Otherwise, we find the ice/regolith boundary:
     74        do ik=1,nsoil-1
     75            if ((icedepth >= layer_PEM(ik)) .and. (icedepth < layer_PEM(ik + 1))) then
     76                iref = ik + 1
     77                exit
     78            endif
     79        enddo
     80        ! And we change the thermal inertia:
     81        do ik = 1,iref - 1
     82            newtherm_i(ig,ik) = max(inert_h2o_ice,prev_thermi(ig,ik))
     83        enddo
     84        ! Transition (based on the equations of thermal conduction):
     85        newtherm_i(ig,iref) = sqrt((layer_PEM(iref) - layer_PEM(iref-1))/(((icedepth - layer_PEM(iref - 1))/newtherm_i(ig,iref - 1)**2) &
     86                              + ((layer_PEM(iref) - icedepth)/inertiedat_PEM(ig,ik)**2) ) )
     87        ! Underlying regolith:
     88        do ik = iref + 1,nsoil
     89            newtherm_i(ig,ik) = inertiedat_PEM(ig,ik)
     90        enddo
     91    endif ! icedepth
     92enddo ! ig
     93
     94return
     95
     96END SUBROUTINE soil_TIfeedback_PEM
     97
     98END MODULE soil_TIfeedback_PEM_mod

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

    r3075 r3076  
    1       subroutine soil_pem_routine(ngrid,nsoil,firstcall, &
    2                therm_i,timestep,tsurf,tsoil)
     1MODULE soil_pem_compute_mod
    32
     3implicit none
    44
    5       use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
    6                           mthermdiff_PEM, thermdiff_PEM, coefq_PEM, &
    7                           coefd_PEM, mu_PEM,alph_PEM,beta_PEM,fluxgeo
    8       use comsoil_h,only: volcapa
    9       implicit none
     5!=======================================================================
     6contains
     7!=======================================================================
     8
     9SUBROUTINE soil_pem_compute(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil)
     10
     11use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
     12use comsoil_h,     only: volcapa
     13
     14implicit none
    1015
    1116!-----------------------------------------------------------------------
    … …  
    2328!  ---------
    2429!  inputs:
    25       integer,intent(in) :: ngrid       ! number of (horizontal) grid-points
    26       integer,intent(in) :: nsoil       ! number of soil layers 
    27       logical,intent(in) :: firstcall    ! identifier for initialization call
    28       real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia [SI]
    29       real,intent(in) :: timestep       ! time step [s]
    30       real,intent(in) :: tsurf(ngrid)   ! surface temperature [K]
     30integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points
     31integer,                      intent(in) :: nsoil     ! number of soil layers 
     32logical,                      intent(in) :: firstcall ! identifier for initialization call
     33real, dimension(ngrid,nsoil), intent(in) :: therm_i  ! thermal inertia [SI]
     34real,                         intent(in) :: timestep  ! time step [s]
     35real, dimension(ngrid),       intent(in) :: tsurf     ! surface temperature [K]
    3136 
    3237! outputs:
    33       real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature [K]
     38real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
    3439! local variables:
    35       integer ig,ik   
     40integer :: ig, ik   
    3641
    3742! 0. Initialisations and preprocessing step
    3843 if (firstcall) then
    39    
    4044! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
    41       do ig=1,ngrid
    42         do ik=0,nsoil-1
    43           mthermdiff_PEM(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa   
    44         enddo
    45       enddo
     45    do ig = 1,ngrid
     46        do ik = 0,nsoil - 1
     47            mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa   
     48        enddo
     49    enddo
    4650
    4751! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
    48       do ig=1,ngrid
    49         do ik=1,nsoil-1
    50       thermdiff_PEM(ig,ik)=((layer_PEM(ik)-mlayer_PEM(ik-1))*mthermdiff_PEM(ig,ik) &
    51                      +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ig,ik-1))  &
    52                          /(mlayer_PEM(ik)-mlayer_PEM(ik-1))
    53         enddo
    54       enddo
     52    do ig = 1,ngrid
     53        do ik = 1,nsoil - 1
     54            thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
     55                                   + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
     56        enddo
     57    enddo
    5558
    5659! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
    57       ! mu_PEM
    58       mu_PEM=mlayer_PEM(0)/(mlayer_PEM(1)-mlayer_PEM(0))
     60    ! mu_PEM
     61    mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
    5962
    60       ! q_{1/2}
    61       coefq_PEM(0)=volcapa*layer_PEM(1)/timestep
    62         ! q_{k+1/2}
    63         do ik=1,nsoil-1
    64           coefq_PEM(ik)=volcapa*(layer_PEM(ik+1)-layer_PEM(ik))                  &
    65                       /timestep
    66         enddo
     63    ! q_{1/2}
     64    coefq_PEM(0) = volcapa*layer_PEM(1)/timestep
     65    ! q_{k+1/2}
     66    do ik = 1,nsoil - 1
     67        coefq_PEM(ik) = volcapa*(layer_PEM(ik + 1) - layer_PEM(ik))/timestep
     68    enddo
    6769
    68       do ig=1,ngrid
    69         ! d_k
    70         do ik=1,nsoil-1
    71           coefd_PEM(ig,ik)=thermdiff_PEM(ig,ik)/(mlayer_PEM(ik)-mlayer_PEM(ik-1))
    72         enddo
    73        
    74         ! alph_PEM_{N-1}
    75         alph_PEM(ig,nsoil-1)=coefd_PEM(ig,nsoil-1)/                      &
    76                        (coefq_PEM(nsoil-1)+coefd_PEM(ig,nsoil-1))
     70    do ig = 1,ngrid
     71        ! d_k
     72        do ik = 1,nsoil - 1
     73            coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik)-mlayer_PEM(ik - 1))
     74        enddo
     75
     76        ! alph_PEM_{N-1}
     77        alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil-1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
    7778        ! alph_PEM_k
    78         do ik=nsoil-2,1,-1
    79           alph_PEM(ig,ik)=coefd_PEM(ig,ik)/(coefq_PEM(ik)+coefd_PEM(ig,ik+1)*    &
    80                                    (1.-alph_PEM(ig,ik+1))+coefd_PEM(ig,ik))
    81         enddo
    82 
     79        do ik = nsoil - 2,1,-1
     80            alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
     81        enddo
    8382      enddo ! of do ig=1,ngrid
    84 
    8583endif ! of if (firstcall)
    8684
    87       IF (.not.firstcall) THEN
    88 !  2. Compute soil temperatures
     85IF (.not. firstcall) THEN
     86! 2. Compute soil temperatures
    8987! First layer:
    90         do ig=1,ngrid
    91           tsoil(ig,1)=(tsurf(ig)+mu_PEM*beta_PEM(ig,1)* & 
    92                                   thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
    93                    (1.+mu_PEM*(1.0-alph_PEM(ig,1))*&
    94                     thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
     88    do ig = 1,ngrid
     89        tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
     90                      (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
    9591
    9692! Other layers:
    97           do ik=1,nsoil-1
    98             tsoil(ig,ik+1)=alph_PEM(ig,ik)*tsoil(ig,ik)+beta_PEM(ig,ik)
    99           enddo
     93        do ik = 1,nsoil - 1
     94                tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
    10095        enddo
    101       ENDIF
     96    enddo
     97endif
    10298
    10399!  2. Compute beta_PEM coefficients (preprocessing for next time step)
    104100! Bottom layer, beta_PEM_{N-1}
    105       do ig=1,ngrid
    106         beta_PEM(ig,nsoil-1)=coefq_PEM(nsoil-1)*tsoil(ig,nsoil)          &
    107                         /(coefq_PEM(nsoil-1)+coefd_PEM(ig,nsoil-1)) &
    108                  +  fluxgeo/(coefq_PEM(nsoil-1)+coefd_PEM(ig,nsoil-1))
    109       enddo
     101do ig = 1,ngrid
     102    beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
     103                             + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
     104enddo
    110105! Other layers
    111       do ik=nsoil-2,1,-1
    112         do ig=1,ngrid
    113           beta_PEM(ig,ik)=(coefq_PEM(ik)*tsoil(ig,ik+1)+                 &
    114                       coefd_PEM(ig,ik+1)*beta_PEM(ig,ik+1))/             &
    115                       (coefq_PEM(ik)+coefd_PEM(ig,ik+1)*(1.0-alph_PEM(ig,ik+1)) &
    116                        +coefd_PEM(ig,ik))
    117         enddo
    118       enddo
     106do ik = nsoil-2,1,-1
     107    do ig = 1,ngrid
     108        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik + 1)*beta_PEM(ig,ik + 1))/ &
     109                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
     110    enddo
     111enddo
    119112
    120       end
     113END SUBROUTINE soil_pem_compute
    121114
     115END MODULE soil_pem_compute_mod

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

    r3075 r3076  
    1       subroutine soil_pem_ini(ngrid,nsoil,therm_i,tsurf,tsoil)
     1MODULE soil_pem_ini_mod
    22
     3implicit none
    34
    4       use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
    5                           mthermdiff_PEM, thermdiff_PEM, coefq_PEM, &
    6                           coefd_PEM, mu_PEM,alph_PEM,beta_PEM,fluxgeo
    7       use comsoil_h,only: volcapa
    8       implicit none
     5!=======================================================================
     6contains
     7!=======================================================================
     8
     9SUBROUTINE soil_pem_ini(ngrid,nsoil,therm_i,tsurf,tsoil)
     10
     11use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
     12use comsoil_h,     only: volcapa
     13
     14implicit none
    915
    1016!-----------------------------------------------------------------------
    1117!  Author: LL
    1218!  Purpose: Compute soil temperature using an implict 1st order scheme, stationnary solution
    13 ! 
    14 !  Note: depths of layers and mid-layers, soil thermal inertia and 
     19!
     20!  Note: depths of layers and mid-layers, soil thermal inertia and
    1521!        heat capacity are commons in comsoil_PEM.h
    1622!-----------------------------------------------------------------------
    … …  
    2228!  ---------
    2329!  inputs:
    24       integer,intent(in) :: ngrid       ! number of (horizontal) grid-points
    25       integer,intent(in) :: nsoil       ! number of soil layers 
    26       real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia [SI]
    27       real,intent(in) :: tsurf(ngrid)   ! surface temperature [K]
     30integer,                      intent(in) :: ngrid   ! number of (horizontal) grid-points
     31integer,                      intent(in) :: nsoil   ! number of soil layers
     32real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI]
     33real, dimension(ngrid),       intent(in) :: tsurf   ! surface temperature [K]
    2834
    2935! outputs:
    30       real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature [K]
     36real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
    3137! local variables:
    32       integer ig,ik,iloop     
     38integer :: ig, ik, iloop
    3339
    3440! 0. Initialisations and preprocessing step
    35    
    3641! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
    37       do ig=1,ngrid
    38         do ik=0,nsoil-1
    39           mthermdiff_PEM(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa   
    40         enddo
    41       enddo
     42do ig = 1,ngrid
     43    do ik = 0,nsoil - 1
     44        mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa
     45    enddo
     46enddo
    4247
    4348! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
    44       do ig=1,ngrid
    45         do ik=1,nsoil-1
    46       thermdiff_PEM(ig,ik)=((layer_PEM(ik)-mlayer_PEM(ik-1))*mthermdiff_PEM(ig,ik) &
    47                      +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ig,ik-1))  &
    48                          /(mlayer_PEM(ik)-mlayer_PEM(ik-1))
    49         enddo
    50       enddo
     49do ig = 1,ngrid
     50    do ik = 1,nsoil - 1
     51        thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
     52                             + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
     53    enddo
     54enddo
    5155
    5256! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
    53       ! mu_PEM
    54       mu_PEM=mlayer_PEM(0)/(mlayer_PEM(1)-mlayer_PEM(0))
     57! mu_PEM
     58mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
    5559
    56       ! q_{1/2}
    57       coefq_PEM(:) = 0.
    58         ! q_{k+1/2}
     60! q_{1/2}
     61coefq_PEM(:) = 0.
     62! q_{k+1/2}
    5963
    60       do ig=1,ngrid
    61         ! d_k
    62         do ik=1,nsoil-1
    63           coefd_PEM(ig,ik)=thermdiff_PEM(ig,ik)/(mlayer_PEM(ik)-mlayer_PEM(ik-1))
    64         enddo
    65        
    66         ! alph_PEM_{N-1}
    67         alph_PEM(ig,nsoil-1)=coefd_PEM(ig,nsoil-1)/                      &
    68                        (coefq_PEM(nsoil-1)+coefd_PEM(ig,nsoil-1))
    69         ! alph_PEM_k
    70         do ik=nsoil-2,1,-1
    71           alph_PEM(ig,ik)=coefd_PEM(ig,ik)/(coefq_PEM(ik)+coefd_PEM(ig,ik+1)*    &
    72                                    (1.-alph_PEM(ig,ik+1))+coefd_PEM(ig,ik))
    73         enddo
     64do ig = 1,ngrid
     65    ! d_k
     66    do ik = 1,nsoil-1
     67        coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
     68    enddo
    7469
    75       enddo ! of do ig=1,ngrid
     70    ! alph_PEM_{N-1}
     71    alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil - 1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
     72    ! alph_PEM_k
     73    do ik = nsoil - 2,1,-1
     74        alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
     75    enddo
     76enddo ! of do ig=1,ngrid
    7677
    77 
    78 
    79 !  1. Compute beta_PEM coefficients
     78!  1. Compute beta_PEM coefficients
    8079! Bottom layer, beta_PEM_{N-1}
    81       do ig=1,ngrid
    82         beta_PEM(ig,nsoil-1)=coefq_PEM(nsoil-1)*tsoil(ig,nsoil)          &
    83                         /(coefq_PEM(nsoil-1)+coefd_PEM(ig,nsoil-1)) &
    84                  +  fluxgeo/(coefq_PEM(nsoil-1)+coefd_PEM(ig,nsoil-1))
    85       enddo
     80do ig = 1,ngrid
     81        beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
     82                                 + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
     83enddo
    8684! Other layers
    87       do ik=nsoil-2,1,-1
    88         do ig=1,ngrid
    89           beta_PEM(ig,ik)=(coefq_PEM(ik)*tsoil(ig,ik+1)+                 &
    90                       coefd_PEM(ig,ik+1)*beta_PEM(ig,ik+1))/             &
    91                       (coefq_PEM(ik)+coefd_PEM(ig,ik+1)*(1.0-alph_PEM(ig,ik+1)) &
    92                        +coefd_PEM(ig,ik))
    93         enddo
    94       enddo
     85do ik = nsoil - 2,1,-1
     86    do ig = 1,ngrid
     87        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik+1)*beta_PEM(ig,ik + 1))/ &
     88                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
     89    enddo
     90enddo
    9591
    9692!  2. Compute soil temperatures
    9793do iloop = 1,10 !just convergence
    98 ! First layer:
    99         do ig=1,ngrid
    100           tsoil(ig,1)=(tsurf(ig)+mu_PEM*beta_PEM(ig,1)* & 
    101                                   thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
    102                    (1.+mu_PEM*(1.0-alph_PEM(ig,1))*&
    103                     thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
    104 ! Other layers:
    105           do ik=1,nsoil-1
    106             tsoil(ig,ik+1)=alph_PEM(ig,ik)*tsoil(ig,ik)+beta_PEM(ig,ik)
    107           enddo
    108         enddo
    109      
    110 enddo ! iloop
    111       end
     94    ! First layer:
     95    do ig = 1,ngrid
     96        tsoil(ig,1)=(tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
     97                   (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
     98    ! Other layers:
     99    do ik = 1,nsoil - 1
     100        tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
     101    enddo
     102enddo
    112103
     104enddo ! iloop
     105
     106END SUBROUTINE soil_pem_ini
     107
     108END MODULE soil_pem_ini_mod

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

    r3075 r3076  
    1       subroutine soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_GCM,
    2      &   TI_GCM,TI_PEM)
     1MODULE soil_settings_PEM_mod
    32
     3implicit none
     4
     5!=======================================================================
     6contains
     7!=======================================================================
     8
     9SUBROUTINE soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_GCM,TI_PEM)
    410
    511!      use netcdf
    6       use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,
    7      &   depth_breccia,depth_bedrock,index_breccia,index_bedrock 
    8       use comsoil_h, only: inertiedat,layer,mlayer, volcapa                       
     12use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, depth_breccia, depth_bedrock, index_breccia, index_bedrock
     13use comsoil_h,     only: inertiedat, layer, mlayer, volcapa
    914
     15use iostart, only: inquire_field_ndims, get_var, get_field, inquire_field, inquire_dimension_length
    1016
    11       use iostart, only: inquire_field_ndims, get_var, get_field,
    12      &                   inquire_field, inquire_dimension_length
     17implicit none
    1318
    14       implicit none
    15 
    16 !======================================================================
     19!=======================================================================
    1720!  Author: LL, based on work by  Ehouarn Millour (07/2006)
    1821!
    1922!  Purpose: Read and/or initialise soil depths and properties
    2023!
    21 ! Modifications: Aug.2010 EM : use NetCDF90 to load variables (enables using
    22 !                      r4 or r8 restarts independently of having compiled
    23 !                      the GCM in r4 or r8)
    24 !                June 2013 TN : Possibility to read files with a time axis
    25 !
     24! Modifications: Aug.2010 EM: use NetCDF90 to load variables (enables using
     25!                              r4 or r8 restarts independently of having compiled
     26!                              the GCM in r4 or r8)
     27!                June 2013 TN: Possibility to read files with a time axis
    2628!
    2729!  The various actions and variable read/initialized are:
    2830!  1. Read/build layer (and midlayer) depth
    29 !  2. Interpolate thermal inertia and temperature on the new grid, if
    30 !     necessary
    31 !======================================================================
     31!  2. Interpolate thermal inertia and temperature on the new grid, if necessary
     32!=======================================================================
    3233
    33 !======================================================================
     34!=======================================================================
    3435!  arguments
    3536!  ---------
    3637!  inputs:
    37       integer,intent(in) :: ngrid       ! # of horizontal grid points
    38       integer,intent(in) :: nslope      ! # of subslope wihtin the mesh
    39       integer,intent(in) :: nsoil_PEM   ! # of soil layers in the PEM
    40       integer,intent(in) :: nsoil_GCM   ! # of soil layers in the GCM
    41       real,intent(in) :: TI_GCM(ngrid,nsoil_GCM,nslope) !  Thermal inertia  in the GCM [SI]
    42       real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope)      ! Thermal inertia   in the PEM [SI]
     38integer,                                 intent(in) :: ngrid     ! # of horizontal grid points
     39integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
     40integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
     41integer,                                 intent(in) :: nsoil_GCM ! # of soil layers in the GCM
     42real, dimension(ngrid,nsoil_GCM,nslope), intent(in) :: TI_GCM    ! Thermal inertia  in the GCM [SI]
    4343
    44 !======================================================================
     44real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Thermal inertia   in the PEM [SI]
     45!=======================================================================
    4546! local variables:
    46       integer ig,iloop,islope   ! loop counters
    47       logical found
    48 
    49       real alpha,lay1 ! coefficients for building layers
    50       real xmin,xmax ! to display min and max of a field
    51      
    52 !======================================================================
    53 
     47integer :: ig, iloop, islope ! loop counters
     48logical :: found
     49real    :: alpha,lay1 ! coefficients for building layers
     50real    :: xmin, xmax ! to display min and max of a field
     51!=======================================================================
    5452! 1. Depth coordinate
    5553! -------------------
    56 
    5754! 1.4 Build mlayer(), if necessary
    58       ! default mlayer distribution, following a power law:
    59       !  mlayer(k)=lay1*alpha**(k-1/2)
    60         lay1=2.e-4
    61         alpha=2
    62         do iloop=0,nsoil_PEM-1
    63           mlayer_PEM(iloop)=lay1*(alpha**(iloop-0.5))
    64         enddo
     55! default mlayer distribution, following a power law:
     56!  mlayer(k)=lay1*alpha**(k-1/2)
     57lay1 = 2.e-4
     58alpha = 2
     59do iloop = 0,nsoil_PEM - 1
     60    mlayer_PEM(iloop) = lay1*(alpha**(iloop - 0.5))
     61enddo
    6562
    6663! 1.5 Build layer(); following the same law as mlayer()
    67       ! Assuming layer distribution follows mid-layer law:
    68       ! layer(k)=lay1*alpha**(k-1)
    69       lay1=sqrt(mlayer_PEM(0)*mlayer_PEM(1))
    70       alpha=mlayer_PEM(1)/mlayer_PEM(0)
    71       do iloop=1,nsoil_PEM
    72         layer_PEM(iloop)=lay1*(alpha**(iloop-1))
    73       enddo
     64! Assuming layer distribution follows mid-layer law:
     65! layer(k)=lay1*alpha**(k-1)
     66lay1 = sqrt(mlayer_PEM(0)*mlayer_PEM(1))
     67alpha = mlayer_PEM(1)/mlayer_PEM(0)
     68do iloop = 1,nsoil_PEM
     69    layer_PEM(iloop) = lay1*(alpha**(iloop - 1))
     70enddo
    7471
    7572! 2. Thermal inertia (note: it is declared in comsoil_h)
    7673! ------------------
    77 
    78       do ig = 1,ngrid
    79         do islope = 1,nslope
    80           do iloop = 1,nsoil_GCM
     74do ig = 1,ngrid
     75    do islope = 1,nslope
     76        do iloop = 1,nsoil_GCM
    8177            TI_PEM(ig,iloop,islope) = TI_GCM(ig,iloop,islope)
    82           enddo
    83           if(nsoil_PEM.gt.nsoil_GCM) then
    84            do iloop = nsoil_GCM+1,nsoil_PEM
    85              TI_PEM(ig,iloop,islope) = TI_GCM(ig,nsoil_GCM,islope)
    86            enddo
    87           endif
    8878        enddo
    89       enddo
    90 
     79        if (nsoil_PEM > nsoil_GCM) then
     80            do iloop = nsoil_GCM + 1,nsoil_PEM
     81                TI_PEM(ig,iloop,islope) = TI_GCM(ig,nsoil_GCM,islope)
     82            enddo
     83        endif
     84    enddo
     85enddo
    9186
    9287! 3. Index for subsurface layering
    9388! ------------------
    94       index_breccia= 1
    95       do iloop = 1,nsoil_PEM-1
    96         if (depth_breccia.ge.layer_PEM(iloop)) then
    97           index_breccia=iloop
    98         else
    99          exit
    100         endif
    101       enddo
     89index_breccia = 1
     90do iloop = 1,nsoil_PEM - 1
     91    if (depth_breccia >= layer_PEM(iloop)) then
     92        index_breccia = iloop
     93    else
     94        exit
     95    endif
     96enddo
    10297
    103       index_bedrock= 1
    104       do iloop = 1,nsoil_PEM-1
    105         if (depth_bedrock.ge.layer_PEM(iloop)) then
    106           index_bedrock=iloop
    107         else
    108          exit
    109         endif
    110       enddo
     98index_bedrock = 1
     99do iloop = 1,nsoil_PEM - 1
     100    if (depth_bedrock >= layer_PEM(iloop)) then
     101        index_bedrock = iloop
     102    else
     103        exit
     104    endif
     105enddo
    111106
     107END SUBROUTINE soil_settings_PEM
    112108
    113       end subroutine soil_settings_PEM
     109END MODULE soil_settings_PEM_mod

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

    r3039 r3076  
    11MODULE time_evol_mod
    22
    3 IMPLICIT NONE 
     3implicit none
    44
    5   integer :: year_bp_ini         !    Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
    6   integer :: dt_pem              !    Time step used by the PEM in Planetary years
    7   real    :: convert_years       !    Conversion ratio from Planetary years to Earth years
    8   real    :: water_ice_criterion !     Percentage of change of the surface of water ice sublimating before stopping the PEM
    9   real    :: co2_ice_criterion   !     Percentage of change of the surface of co2 ice sublimating before stopping the PEM
    10   real    :: ps_criterion        !    Percentage of change of averaged surface pressure before stopping the PEM
    11   integer :: Max_iter_pem        !    Maximal number of iteration when converging to a steady state, read in evol.def
    12   logical :: evol_orbit_pem      !    True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
    13   logical :: var_obl             !     True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
    14   logical :: var_ecc             !    True if we want the PEM to follow obl_ecc_lsp.asc parameters for excenticity
    15   logical :: var_lsp             !    True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
     5integer :: year_bp_ini         ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
     6integer :: dt_pem              ! Time step used by the PEM in Planetary years
     7real    :: convert_years       ! Conversion ratio from Planetary years to Earth years
     8real    :: water_ice_criterion ! Percentage of change of the surface of water ice sublimating before stopping the PEM
     9real    :: co2_ice_criterion   ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
     10real    :: ps_criterion        ! Percentage of change of averaged surface pressure before stopping the PEM
     11integer :: Max_iter_pem        ! Maximal number of iteration when converging to a steady state, read in evol.def
     12logical :: evol_orbit_pem      ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
     13logical :: var_obl             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
     14logical :: var_ecc             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for excenticity
     15logical :: var_lsp             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
    1616
    1717END MODULE time_evol_mod

  • trunk/LMDZ.MARS/deftank/pem/launch_pem.sh

    r3068 r3076  
    168168    mv startfi.nc startfi_evol.nc
    169169    if [ -f "start.nc" ]; then
    170         cp start.nc start_evol.nc
     170        mv start.nc start_evol.nc
    171171    elif [ -f "start1D.txt" ]; then
    172         cp start1D.txt start1D_evol.txt
     172        mv start1D.txt start1D_evol.txt
    173173    fi
    174174    ./$exePEM > out_runPEM${iPEM} 2>&1
Note: See TracChangeset for help on using the changeset viewer.