Ignore:
Timestamp:
Jun 12, 2023, 4:38:28 PM (19 months ago)
Author:
romain.vande
Message:

Mars PEM :

Adapt PEM to 1d runs.
Cleaning of names and unused variables.
Correct minor errors.
Adapt and correct reshape_xios_output utilitary for 1d diagfi output.

RV

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
3 added
5 deleted
16 edited

Legend:

Unmodified
Added
Removed

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

    r2888 r2980  
    11!
    2 ! $Id: abort_pem.F 1425 2010-09-02 13:45:23Z lguez $
     2! $Id: abort_pem.F
    33!
    44c
    … …  
    4444#endif
    4545      call getin_dump
    46 c     call histclo(2)
    47 c     call histclo(3)
    48 c     call histclo(4)
    49 c     call histclo(5)
    5046      write(lunout,*) 'Stopping in ', modname
    5147      write(lunout,*) 'Reason = ',message

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

    r2944 r2980  
    1313!!!
    1414!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    15 
    16 
    1715  subroutine ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
    1816
    … …  
    2523  end subroutine ini_adsorption_h_PEM
    2624
     25!!! -----------------------------------------------
     26
    2727  subroutine end_adsorption_h_PEM
    2828
    … …  
    3232  end subroutine end_adsorption_h_PEM
    3333
    34 
    35 
    36 
    37 
     34!!! -----------------------------------------------
    3835
    3936  subroutine regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
    … …  
    354351                                   theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o)
    355352
    356 
    357 
    358353! 2.  we compute the mass of co2 adsorded in each layer of the meshes 
    359354

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

    r2945 r2980  
    22
    33implicit none
    4   integer, parameter :: nsoilmx_PEM = 27               ! number of layers in the PEM
    5   real,save,allocatable,dimension(:) :: layer_PEM      ! soil layer depths   [m]
    6   real,save,allocatable,dimension(:) :: mlayer_PEM     ! soil mid-layer depths [m]
    7   real,save,allocatable,dimension(:,:,:) :: TI_PEM     ! soil thermal inertia [SI]
     4  integer, parameter :: nsoilmx_PEM = 27                 ! number of layers in the PEM
     5  real,save,allocatable,dimension(:) :: layer_PEM        ! soil layer depths   [m]
     6  real,save,allocatable,dimension(:) :: mlayer_PEM       ! soil mid-layer depths [m]
     7  real,save,allocatable,dimension(:,:,:) :: TI_PEM       ! soil thermal inertia [SI]
    88  real,save,allocatable,dimension(:,:) :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI]
    99  ! variables (FC: built in firstcall in soil.F)
    10   REAL,SAVE,ALLOCATABLE :: tsoil_PEM(:,:,:)       ! sub-surface temperatures [K]
    11   real,save,allocatable :: mthermdiff_PEM(:,:)  ! (FC) mid-layer thermal diffusivity [SI]
    12   real,save,allocatable :: thermdiff_PEM(:,:)   ! (FC) inter-layer thermal diffusivity [SI]
    13   real,save,allocatable :: coefq_PEM(:)         ! (FC) q_{k+1/2} coefficients [SI]
    14   real,save,allocatable :: coefd_PEM(:,:)       ! (FC) d_k coefficients [SI]
    15   real,save,allocatable :: alph_PEM(:,:)      ! (FC) alpha_k coefficients [SI]
    16   real,save,allocatable :: beta_PEM(:,:)      ! beta_k coefficients [SI]
    17   real,save :: mu_PEM                           ! mu coefficient [SI]
    18   real,save :: fluxgeo                               ! Geothermal flux [W/m^2]
    19   real,save :: depth_breccia                         ! Depth at which we have breccia [m]
    20   real,save :: depth_bedrock                         ! Depth at which we have bedrock [m]
    21   integer,save :: index_breccia                         ! last index of the depth grid before having breccia
    22   integer,save :: index_bedrock                         ! last index of the depth grid before having bedrock
     10  REAL,SAVE,ALLOCATABLE :: tsoil_PEM(:,:,:)              ! sub-surface temperatures [K]
     11  real,save,allocatable :: mthermdiff_PEM(:,:)           ! (FC) mid-layer thermal diffusivity [SI]
     12  real,save,allocatable :: thermdiff_PEM(:,:)            ! (FC) inter-layer thermal diffusivity [SI]
     13  real,save,allocatable :: coefq_PEM(:)                  ! (FC) q_{k+1/2} coefficients [SI]
     14  real,save,allocatable :: coefd_PEM(:,:)                ! (FC) d_k coefficients [SI]
     15  real,save,allocatable :: alph_PEM(:,:)                 ! (FC) alpha_k coefficients [SI]
     16  real,save,allocatable :: beta_PEM(:,:)                 ! beta_k coefficients [SI]
     17  real,save :: mu_PEM                                    ! mu coefficient [SI]
     18  real,save :: fluxgeo                                   ! Geothermal flux [W/m^2]
     19  real,save :: depth_breccia                             ! Depth at which we have breccia [m]
     20  real,save :: depth_bedrock                             ! Depth at which we have bedrock [m]
     21  integer,save :: index_breccia                          ! last index of the depth grid before having breccia
     22  integer,save :: index_bedrock                          ! last index of the depth grid before having bedrock
    2323  LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
    24   real,save,allocatable,dimension(:) :: water_reservoir      ! Large reserve of water   [kg/m^2]
    25   real,save :: water_reservoir_nom
    26   logical, save :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
     24  real,save,allocatable,dimension(:) :: water_reservoir  ! Large reserve of water   [kg/m^2]
     25  real,save :: water_reservoir_nom                       ! Nominal value for the large reserve of water   [kg/m^2]
     26  logical, save :: reg_thprop_dependp                    ! thermal properites of the regolith vary with the surface pressure
    2727
    2828contains

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

    r2963 r2980  
    11MODULE conf_pem_mod
     2
     3!=======================================================================
     4!
     5! Purpose: Read the parameter for a PEM run in run_pem.def
     6!
     7! Author: RV
     8!=======================================================================
    29
    310IMPLICIT NONE

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

    r2893 r2980  
    3030
    3131!   INPUT
    32   INTEGER, intent(in) :: ngrid                  ! # of grid physical grid points
    33   REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
    34   REAL,    intent(in) :: qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
    35   REAL,    intent(in) :: ini_surf
    36   REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
     32  INTEGER, intent(in) :: ngrid                         ! # of grid physical grid points
     33  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
     34  REAL,    intent(in) :: qsurf(ngrid,nslope)           ! physical point field : Actual density of water ice
     35  REAL,    intent(in) :: ini_surf                      ! Initial surface of h2o ice that was sublimating
     36  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope) ! Grid point that initialy were covered by h2o_ice
    3737
    3838!   OUTPUT
    … …  
    4141!   local:
    4242!   -----
    43   INTEGER :: i,islope                    ! Loop
     43  INTEGER :: i,islope   ! Loop
    4444  REAL :: present_surf  ! Initial/Actual surface of water ice
    4545
    … …  
    7777! ------------------------------------------------------------------------------------------------
    7878
    79 
    80 SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_h2o_ice,global_ave_press_GCM,global_ave_press_new,nslope)
     79SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_co2_ice,global_ave_press_GCM,global_ave_press_new,nslope)
    8180
    8281  USE temps_mod_evol, ONLY: co2_ice_criterion,ps_criterion
    … …  
    9897  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
    9998  REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
    100   REAL,    intent(in) :: ini_surf
    101   REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
    102   REAL,    intent(in) :: global_ave_press_GCM
    103   REAL,    intent(in) :: global_ave_press_new
     99  REAL,    intent(in) :: ini_surf                      ! Initial surface of co2 ice that was sublimating
     100  REAL,    intent(in) :: initial_co2_ice(ngrid,nslope) ! Grid point that initialy were covered by co2_ice
     101  REAL,    intent(in) :: global_ave_press_GCM          ! Planet average pressure from the GCM start files
     102  REAL,    intent(in) :: global_ave_press_new          ! Planet average pressure from the PEM computations
    104103
    105104!   OUTPUT
    … …  
    120119  do i=1,ngrid
    121120   do islope=1,nslope
    122       if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
     121      if (initial_co2_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
    123122         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
    124123      endif

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

    r2961 r2980  
    3737  end subroutine end_ice_table_porefilling
    3838
    39 
    40 
     39!!! --------------------------------------
    4140
    4241   SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
    … …  
    134133      enddo
    135134
    136 !=======================================================================
    137135      RETURN
    138136#endif
    139137      END
    140138
     139!!! --------------------------------------
    141140
    142141   SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
    … …  
    189188real :: pvap2rho = 18.e-3/8.314
    190189!
    191 
    192 
    193 
    194 
    195190 
    196191! 0. Compute constriction over the layer
    … …  
    300295end subroutine
    301296
    302 
    303 
     297!!! --------------------------------------
    304298
    305299SUBROUTINE constriction(porefill,eta)
    … …  
    318312        !!!
    319313
    320 
    321314        if (porefill.le.0.) eta = 1.
    322315        if ((porefill.gt.0.) .and.(porefill.lt.1.)) then
    … …  
    327320        end subroutine
    328321
    329 
    330322end module

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

    r2888 r2980  
    11SUBROUTINE info_run_PEM(year_iter, criterion_stop)
     2
     3!=======================================================================
     4!
     5! Purpose: Write in a file called info_run_PEM.txt the reason why the PEM did stop and the number of extrapolation year done
     6!
     7! Author: RV
     8!=======================================================================
    29
    310   IMPLICIT NONE
    411
    5   INTEGER, intent(in) :: year_iter, criterion_stop                  ! # of year and reason to stop
     12  INTEGER, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
    613
    714  logical :: exist
    8 
    915
    1016  inquire(file='info_run_PEM.txt', exist=exist)

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

    r2855 r2980  
    1111
    1212  real,save,allocatable :: yearlask(:)          ! year before present from Laskar et al. Tab
    13   real,save,allocatable :: oblask(:)            ! obliquity [deg]
    14   real,save,allocatable :: exlask(:)            ! excentricity[deg]
     13  real,save,allocatable :: oblask(:)            ! obliquity    [deg]
     14  real,save,allocatable :: exlask(:)            ! excentricity [deg]
    1515  real,save,allocatable :: lsplask(:)           ! ls perihelie [deg]
    1616  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

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

    r2855 r2980  
    4444          write(*,*)"read_data_GCM: Le champ <time_counter> est absent"
    4545          write(*,*)"read_data_GCM: J essaie <Time>"
     46          ierr = nf90_inq_varid (fID, "Time", vID)
    4647          IF (ierr .NE. nf90_noerr) THEN
    4748            write(*,*)"read_data_GCM: Le champ <Time> est absent"
    4849            write(*,*)trim(nf90_strerror(ierr))
    49             CALL ABORT_gcm("read_data_GCM", "", 1)
     50            CALL ABORT_gcm("nb_time_step_GCM", "", 1)
    5051          ENDIF
    5152          ! Get the length of the "Time" dimension
    5253          ierr = nf90_inq_dimid(fID,"Time",vID)
    5354          ierr = nf90_inquire_dimension(fID,vID,len=timelen)         
    54         ENDIF
     55        ELSE
    5556        ! Get the length of the "time_counter" dimension
    5657        ierr = nf90_inq_dimid(fID,"time_counter",vID)
    5758        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
     59        ENDIF
    5860      ELSE   
    5961        ! Get the length of the "temps" dimension

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

    r2963 r2980  
    11      MODULE orbit_param_criterion_mod
     2
     3!=======================================================================
     4!
     5! Purpose: Compute the maximum nimber of iteration the PEM can do based
     6! on the stopping criterion given by the orbital parameters
     7!
     8! Author: RV
     9!=======================================================================
    210
    311      IMPLICIT NONE
    … …  
    3442
    3543      real :: Year                      ! Year of the simulation
    36       real :: timeperi_ls                       ! Year of the simulation
    37       integer nlask,ilask!,last_ilask !Loop variable
    38       parameter (nlask = 20001)
     44      real :: timeperi_ls               ! Year of the simulation
     45      integer nlask,ilask!,last_ilask   ! Loop variable
     46      parameter (nlask = 20001)         ! Number of line in Laskar file
    3947
    4048      real max_change_obl,max_change_ex,max_change_lsp ! Percentage of change that is considered to be acceptible
    4149
    42       real max_obl,max_ex,max_lsp !Maximum value of orbit param given the acceptable percentage
    43       real min_obl,min_ex,min_lsp !Maximum value of orbit param given the acceptable percentage
    44 
    45       real max_obl_iter,max_ex_iter,max_lsp_iter !Maximum year iteration before reaching an unacceptable value
    46 
    47       logical obl_not_found, ex_not_found,lsp_not_found !Loop variable (first call)
    48       logical max_lsp_modulo,min_lsp_modulo
     50      real max_obl,max_ex,max_lsp       ! Maximum value of orbit param given the acceptable percentage
     51      real min_obl,min_ex,min_lsp       ! Maximum value of orbit param given the acceptable percentage
     52
     53      real max_obl_iter,max_ex_iter,max_lsp_iter        ! Maximum year iteration before reaching an unacceptable value
     54
     55      logical obl_not_found, ex_not_found,lsp_not_found ! Loop variable (first call)
     56      logical max_lsp_modulo,min_lsp_modulo             ! Variable to check if we are close to the 360 degree modulo for lsp parameter
     57
     58      real xa,xb,ya,yb,yc
    4959
    5060      ! **********************************************************************
    … …  
    5262      ! **********************************************************************
    5363
    54           Year=year_bp_ini+year_PEM
    55           timeperi_ls=timeperi*360/2/pi
    56 
    57           call ini_lask_param_mod(nlask)
     64          Year=year_bp_ini+year_PEM            ! We compute the current year
     65          timeperi_ls=timeperi*360/2/pi        ! We convert in degree
     66
     67          call ini_lask_param_mod(nlask)       ! Allocation of variables
    5868
    5969          print *, "orbit_param_criterion, Year in pem.def=", year_bp_ini
    … …  
    6373          print *, "orbit_param_criterion, Current ex=", e_elips
    6474          print *, "orbit_param_criterion, Current lsp=", timeperi_ls
     75
     76! We read the file
    6577
    6678          open(73,file='ob_ex_lsp.asc')
    … …  
    139151        max_lsp_iter=999999999999
    140152
     153!--------------------------------
     154
     155        yc=max_obl
     156        yc=min_obl
     157
     158!--------------------------------
     159
    141160        do ilask=last_ilask,1,-1
    142            if((oblask(ilask).GT.max_obl).and. obl_not_found ) then
    143               max_obl_iter=(max_obl-oblask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
    144                                / (oblask(ilask)-oblask(ilask-1))+yearlask(ilask-1)-Year
     161
     162          xa=yearlask(ilask)
     163          xb=yearlask(ilask-1)
     164! Linear interpolation to find the maximum number of iteration
     165           if((oblask(ilask-1).GT.max_obl).and. obl_not_found ) then
     166              ya=oblask(ilask)
     167              yb=oblask(ilask-1)
     168              yc=max_obl
     169              max_obl_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year
    145170              obl_not_found=.FALSE.
    146            elseif((oblask(ilask).LT.min_obl).and. obl_not_found ) then
    147               max_obl_iter=(min_obl-oblask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
    148                                / (oblask(ilask)-oblask(ilask-1))+yearlask(ilask-1)-Year
     171           elseif((oblask(ilask-1).LT.min_obl).and. obl_not_found ) then
     172              ya=oblask(ilask)
     173              yb=oblask(ilask-1)
     174              yc=min_obl
     175              max_obl_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year
    149176              obl_not_found=.FALSE.
    150177           endif
    151            if((exlask(ilask).GT.max_ex).and. ex_not_found ) then
    152               max_ex_iter=(max_ex-exlask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
    153                                / (exlask(ilask)-exlask(ilask-1))+yearlask(ilask-1)-Year
     178           if((exlask(ilask-1).GT.max_ex).and. ex_not_found ) then
     179              ya=exlask(ilask)
     180              yb=exlask(ilask-1)
     181              yc=max_ex
     182              max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year
    154183              ex_not_found=.FALSE.
    155            elseif((exlask(ilask).LT.min_ex ).and. ex_not_found ) then
    156               max_ex_iter=(min_ex-exlask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
    157                                / (exlask(ilask)-exlask(ilask-1))+yearlask(ilask-1)-Year
     184           elseif((exlask(ilask-1).LT.min_ex ).and. ex_not_found ) then
     185              ya=exlask(ilask)
     186              yb=exlask(ilask-1)
     187              yc=min_ex
     188              max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year
    158189              ex_not_found=.FALSE.
    159190           endif
    160191           if(max_lsp_modulo .and. lsplask(ilask).lt.180) lsplask(ilask)=lsplask(ilask)+360.
    161192           if(min_lsp_modulo .and. lsplask(ilask).gt.180) lsplask(ilask)=lsplask(ilask)-360.
    162            if((lsplask(ilask).GT.max_lsp).and. lsp_not_found ) then
    163               max_lsp_iter=(max_lsp-lsplask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
    164                                / (lsplask(ilask)-lsplask(ilask-1))+yearlask(ilask-1)-Year
     193           if((lsplask(ilask-1).GT.max_lsp).and. lsp_not_found ) then
     194              ya=lsplask(ilask)
     195              yb=lsplask(ilask-1)
     196              yc=max_lsp
     197              max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year
    165198              lsp_not_found=.FALSE.
    166            elseif((lsplask(ilask).LT.min_lsp ).and. lsp_not_found ) then
    167               max_lsp_iter=(min_lsp-lsplask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
    168                                / (lsplask(ilask)-lsplask(ilask-1))+yearlask(ilask-1)-Year
     199           elseif((lsplask(ilask-1).LT.min_lsp ).and. lsp_not_found ) then
     200              ya=lsplask(ilask)
     201              yb=lsplask(ilask-1)
     202              yc=min_lsp
     203              max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year
    169204              lsp_not_found=.FALSE.
    170205           endif

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

    r2963 r2980  
    5454                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
    5555                        nf90_inquire_dimension,nf90_close
    56 
     56#ifndef CPP_1D
     57      USE iniphysiq_mod, ONLY: iniphysiq
    5758! For phyredem :
    5859      USE control_mod, ONLY: iphysiq, day_step,nsplit_phys
    59       USE iniphysiq_mod, ONLY: iniphysiq
     60#else
     61      use time_phylmdz_mod, only: daysec, iphysiq,day_step
     62#endif
    6063      USE logic_mod, ONLY: iflag_phys
    6164#ifndef CPP_STD
    … …  
    8588      use constants_marspem_mod,only: alpha_clap_co2,beta_clap_co2, alpha_clap_h2o,beta_clap_h2o, &
    8689                                      m_co2,m_noco2
     90      use evol_co2_ice_s_mod, only: evol_co2_ice_s
     91      use evol_h2o_ice_s_mod, only: evol_h2o_ice_s
    8792
    8893!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
    … …  
    105110      use soil_thermalproperties_mod, only: update_soil_thermalproperties
    106111
     112#ifdef CPP_1D
     113      use regular_lonlat_mod, only: init_regular_lonlat
     114      use physics_distribution_mod, only: init_physics_distribution
     115      use mod_grid_phy_lmdz, only : regular_lonlat
     116      use init_phys_1d_mod, only : init_phys_1d
     117#endif
     118
    107119
    108120  IMPLICIT NONE
    … …  
    114126  PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
    115127
    116   include "comdissnew.h"
    117128  include "comgeom.h"
    118129  include "iniprint.h"
    … …  
    177188      LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
    178189      LOGICAL :: STOPPING_co2     ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
    179       LOGICAL :: STOPPING_1_co2   ! Logical : is there still water ice to sublimate?
    180190      LOGICAL :: STOPPING_pressure
    181191      INTEGER :: criterion_stop  ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
    182192
    183       real,save ::     A , B, mmean        ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
    184       real ,allocatable :: vmr_co2_gcm(:,:) ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
     193      real,save ::     A , B, mmean              ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
     194      real ,allocatable :: vmr_co2_gcm(:,:)      ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
    185195      real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times  co2 volume mixing ratio used in the PEM
    186196      real ,allocatable :: q_co2_PEM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg]
    187       REAL ,allocatable :: q_h2o_PEM_phys(:,:)  ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]
     197      REAL ,allocatable :: q_h2o_PEM_phys(:,:)   ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]
    188198      integer :: timelen                         ! # time samples
    189199      REAL :: ave                                ! intermediate varibale to compute average
    … …  
    194204
    195205!!!!!!!!!!!!!!!!!!!!!!!! SLOPE
    196       REAL ::              watercap_saved                          ! Value saved at the previous time step
    197       REAL , dimension(:,:), allocatable :: min_co2_ice_1                ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
    198       REAL , dimension(:,:), allocatable :: min_co2_ice_2                ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
    199       REAL , dimension(:,:), allocatable :: min_h2o_ice_1                ! ngrid field : minimum of water ice at each point for the first year [kg/m^2]
    200       REAL , dimension(:,:), allocatable :: min_h2o_ice_2                ! ngrid field : minimum of water ice at each point for the second year [kg/m^2]
    201       REAL , dimension(:,:,:), allocatable :: co2_ice_GCM          ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
    202       REAL, dimension(:,:),allocatable  :: initial_co2_ice_sublim    ! physical point field : Logical array indicating sublimating point of co2 ice
    203       REAL, dimension(:,:),allocatable  :: initial_h2o_ice           ! physical point field : Logical array indicating if there is water ice at initial state
    204       REAL, dimension(:,:),allocatable  :: initial_co2_ice           ! physical point field : Logical array indicating if there is co2 ice at initial state
    205       REAL, dimension(:,:),allocatable  :: tendencies_co2_ice              ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
    206       REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_ini          ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
    207       REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice              ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
    208       REAL , dimension(:,:), allocatable :: flag_co2flow(:,:)   !(ngrid,nslope)          ! To flag where there is a glacier flow
    209       REAL , dimension(:), allocatable :: flag_co2flow_mesh(:)  !(ngrid)          ! To flag where there is a glacier flow
     206      REAL                               :: watercap_saved          ! Value saved at the previous time step
     207      REAL, dimension(:,:),  allocatable :: min_co2_ice_1           ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
     208      REAL, dimension(:,:),  allocatable :: min_co2_ice_2           ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
     209      REAL, dimension(:,:),  allocatable :: min_h2o_ice_1           ! ngrid field : minimum of water ice at each point for the first year [kg/m^2]
     210      REAL, dimension(:,:),  allocatable :: min_h2o_ice_2           ! ngrid field : minimum of water ice at each point for the second year [kg/m^2]
     211      REAL, dimension(:,:,:),allocatable :: co2_ice_GCM             ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
     212      REAL, dimension(:,:),  allocatable :: initial_co2_ice_sublim  ! physical point field : Logical array indicating sublimating point of co2 ice
     213      REAL, dimension(:,:),  allocatable :: initial_h2o_ice         ! physical point field : Logical array indicating if there is water ice at initial state
     214      REAL, dimension(:,:),  allocatable :: initial_co2_ice         ! physical point field : Logical array indicating if there is co2 ice at initial state
     215      REAL, dimension(:,:),  allocatable :: tendencies_co2_ice      ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
     216      REAL, dimension(:,:),  allocatable :: tendencies_co2_ice_ini  ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
     217      REAL, dimension(:,:),  allocatable :: tendencies_h2o_ice      ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
     218      REAL, dimension(:,:),  allocatable :: flag_co2flow(:,:)       !(ngrid,nslope): Flag where there is a glacier flow
     219      REAL, dimension(:),    allocatable :: flag_co2flow_mesh(:)    !(ngrid)       : Flag where there is a glacier flow
    210220
    211221!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL
    … …  
    213223     REAL, ALLOCATABLE :: tsurf_ave(:,:)       ! Physic x SLOPE field : Averaged Surface Temperature [K]
    214224     REAL, ALLOCATABLE :: tsoil_ave(:,:,:)   ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
    215      REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:) ! Physics x SLOPE field : Averaged Soil Temperature  during 1st year [K]
    216      REAL, ALLOCATABLE :: tsoil_ave_PEM_yr1(:,:,:)
    217225     REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:)    ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
    218226
    … …  
    225233     REAL,ALLOCATABLE  :: Tsurf_locslope(:)   ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
    226234     REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
    227      REAL,ALLOCATABLE  :: watersurf_density_ave(:,:)          ! Physic  x Slope, water surface density, yearly averaged [kg/m^3]
     235     REAL,ALLOCATABLE  :: watersurf_density_ave(:,:)            ! Physic  x Slope, water surface density, yearly averaged [kg/m^3]
    228236     REAL,ALLOCATABLE  :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
    229237     REAL,ALLOCATABLE  :: watersoil_density_PEM_ave(:,:,:)          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
    230      REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)                     ! Surface temperature saved from previous time step [K]
    231      REAL, ALLOCATABLE :: delta_co2_adsorbded(:)                         ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    232      REAL, ALLOCATABLE :: delta_h2o_adsorbded(:)                         ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
    233      REAL :: totmassco2_adsorbded                                           ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
    234      REAL :: totmassh2o_adsorbded                                           ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
    235      LOGICAL :: bool_sublim                                              ! logical to check if there is sublimation or not
     238     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)                ! Surface temperature saved from previous time step [K]
     239     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)                    ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
     240     REAL, ALLOCATABLE :: delta_h2o_adsorbded(:)                    ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
     241     REAL :: totmassco2_adsorbded                                   ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
     242     REAL :: totmassh2o_adsorbded                                   ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
     243     LOGICAL :: bool_sublim                                         ! logical to check if there is sublimation or not
    236244     
    237245
    … …  
    249257     REAL :: albedo_h2o_frost                        ! albedo of h2o frost
    250258     REAL,ALLOCATABLE :: co2ice(:)                   ! Physics: co2 ice mesh averaged [kg/m^2]
     259#endif
     260
     261#ifdef CPP_1D
     262INTEGER :: nsplit_phys
     263integer,parameter:: jjm_value=jjm-1
     264integer :: ierr
     265#else
     266integer,parameter:: jjm_value=jjm
    251267#endif
    252268
    … …  
    277293      year_day=669
    278294      daysec=88775.
    279       dtphys=0
    280295      timestep=year_day*daysec/year_step
    281296
    … …  
    288303
    289304!----------------------------READ run.def ---------------------
     305
     306#ifndef CPP_1D
     307      dtphys=0
    290308      CALL conf_gcm( 99, .TRUE. )
     309      call infotrac_init
     310      nq=nqtot
     311     allocate(q(ip1jmp1,llm,nqtot))
     312#else
     313      ! load tracer names from file 'traceur.def'
     314        open(90,file='traceur.def',status='old',form='formatted',&
     315            iostat=ierr)
     316        if (ierr.ne.0) then
     317          write(*,*) 'Cannot find required file "traceur.def"'
     318          write(*,*) ' If you want to run with tracers, I need it'
     319          write(*,*) ' ... might as well stop here ...'
     320          stop
     321        else
     322          write(*,*) "pem1d: Reading file traceur.def"
     323          ! read number of tracers:
     324          read(90,*,iostat=ierr) nq
     325          print *, "nq",nq
     326          nqtot=nq ! set value of nqtot (in infotrac module) as nq
     327          if (ierr.ne.0) then
     328            write(*,*) "pem1d: error reading number of tracers"
     329            write(*,*) "   (first line of traceur.def) "
     330            stop
     331          endif
     332          if (nq<1) then
     333            write(*,*) "pem1d: error number of tracers"
     334            write(*,*) "is nq=",nq," but must be >=1!"
     335            stop
     336          endif
     337        endif
     338     nq=nqtot
     339     allocate(q(ip1jmp1,llm,nqtot))
     340     allocate(ap(nlayer+1))
     341     allocate(bp(nlayer+1))
     342     call init_phys_1d(llm,nqtot,vcov,ucov, &
     343                    teta,q,ps, time_0,ap,bp)   
     344     pi=2.E+0*asin(1.E+0)
     345     g=3.72
     346     nsplit_phys=1
     347#endif
    291348      CALL conf_pem     
    292349 
    293       call infotrac_init
    294       nq=nqtot
    295 
    296350!------------------------
    297351
    … …  
    306360!----------------------------READ start.nc ---------------------
    307361
    308      allocate(q(ip1jmp1,llm,nqtot))
     362
     363#ifndef CPP_1D
    309364     CALL dynetat0(FILE_NAME_start,vcov,ucov, &
    310365                    teta,q,masse,ps,phis, time_0)
     366#endif
    311367
    312368     allocate(ps_start_GCM(ngrid))
     369#ifndef CPP_1D
    313370     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
    314371
    315372     CALL iniconst !new
    316373     CALL inigeom
     374
    317375     allocate(ap(nlayer+1))
    318376     allocate(bp(nlayer+1))
    … …  
    324382     status =nf90_close(ncid)
    325383
     384#else
     385
     386     ps_start_GCM(1)=ps(1)
     387
     388#endif
     389
     390
     391#ifndef CPP_1D
    326392     CALL iniphysiq(iim,jjm,llm, &
    327393          (jjm-1)*iim+2,comm_lmdz, &
    … …  
    329395          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
    330396          iflag_phys)
     397#endif
    331398
    332399! In the gcm, these values are given to the physic by the dynamic.
    … …  
    470537     allocate(tsurf_ave_yr1(ngrid,nslope))
    471538     allocate(tsurf_ave(ngrid,nslope))
    472      allocate(tsoil_ave_yr1(ngrid,nsoilmx,nslope))
    473539     allocate(tsoil_ave(ngrid,nsoilmx,nslope))
    474540     allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
    … …  
    480546     allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
    481547
    482      allocate(tsoil_ave_PEM_yr1(ngrid,nsoilmx_PEM,nslope))
    483548     allocate(Tsurfave_before_saved(ngrid,nslope))
    484549     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    … …  
    491556     print *, "Downloading data Y1..."
    492557
    493      call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1,&   
    494                        tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, &     
     558     call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1,&   
     559                       tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, &     
    495560                       watersurf_density_ave,watersoil_density_timeseries)
    496561
    … …  
    500565     print *, "Downloading data Y2"
    501566
    502      call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
     567     call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
    503568                  tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, &     
    504569                  watersurf_density_ave,watersoil_density_timeseries)
    … …  
    529594      call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
    530595      DO l=1,nsoilmx
    531         tsoil_ave_PEM_yr1(:,l,:)=tsoil_ave_yr1(:,l,:)
    532596        tsoil_PEM(:,l,:)=tsoil_ave(:,l,:)
    533597        tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:)
    … …  
    535599      ENDDO
    536600      DO l=nsoilmx+1,nsoilmx_PEM
    537         tsoil_ave_PEM_yr1(:,l,:)=tsoil_ave_yr1(:,nsoilmx,:)
    538601        tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:)
    539602        watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:)
    540603      ENDDO
    541604      watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
    542       deallocate(tsoil_ave_yr1)
    543605      deallocate(tsoil_ave)
    544606      deallocate(tsoil_GCM_timeseries)
    … …  
    656718!---------------------------- Read the PEMstart ---------------------
    657719
    658       call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_PEM_yr1,tsoil_PEM,porefillingice_depth,porefillingice_thickness, &
     720      call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth,porefillingice_thickness, &
    659721      tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,&
    660722      tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,&
    … …  
    687749     write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
    688750     endif ! adsorption
    689      deallocate(tsoil_ave_PEM_yr1) 
    690751     deallocate(tsurf_ave_yr1)
    691752
    … …  
    702763!    I_h Read the PEMstar
    703764!    I_i Compute orbit criterion
     765
    704766#ifndef CPP_STD
    705767     CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
    … …  
    813875      print *, "Evolution of h2o ice"
    814876     
    815      call evol_h2o_ice_s_slope(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
     877     call evol_h2o_ice_s(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope)
    816878     
    817879     DO islope=1, nslope
    … …  
    822884
    823885      print *, "Evolution of co2 ice"
    824      call evol_co2_ice_s_slope(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
     886     call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
    825887
    826888!------------------------
    … …  
    839901                         global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
    840902      endif
    841 
    842 
    843903
    844904     DO islope=1, nslope
    … …  
    10311091        criterion_stop=2
    10321092      endif
    1033       if (STOPPING_1_co2) then
    1034         print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
    1035         criterion_stop=2
    1036       endif
    10371093      if (STOPPING_pressure) then
    10381094        print *, "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
    … …  
    10441100      endif
    10451101
    1046      if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2 .or. STOPPING_pressure)  then
     1102     if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure)  then
    10471103        exit
    10481104      else
    … …  
    11811237       enddo
    11821238     enddo
    1183      
     1239   
    11841240!------------------------
    11851241     if(evol_orbit_pem) then
    … …  
    12001256
    12011257     allocate(p(ip1jmp1,nlayer+1))
     1258#ifndef CPP_1D
    12021259     CALL pression (ip1jmp1,ap,bp,ps,p)
    12031260     CALL massdair(p,masse)
    … …  
    12081265                time_0,vcov,ucov,teta,q,masse,ps)
    12091266      print *, "restart_evol.nc has been written"
     1267
     1268#else
     1269     DO nnq = 1, nqtot
     1270        call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq)
     1271     ENDDO
     1272#endif
    12101273
    12111274! III_b.2 WRITE restartfi.nc

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

    r2961 r2980  
    1    SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, ice_table_thickness, &
     1   SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
    22                      tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
    33                      global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
    … …  
    3535  real,intent(in) :: co2ice(ngrid,nslope)                     ! co2 ice amount [kg/m^2]
    3636  real,intent(in) :: waterice(ngrid,nslope)                   ! water ice amount [kg/m^2]
    37   real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope)   ! soil temperature during 1st year [K]
    3837  real, intent(in) :: watersurf_ave(ngrid,nslope)             ! surface water ice density, yearly averaged  (kg/m^3)
    3938  real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)! surface water ice density, yearly averaged (kg/m^3)

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

    r2963 r2980  
    1 !
    2 ! $Id $
    3 !
    41SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries, &
    52             min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
    … …  
    2825  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
    2926  INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
    30   INTEGER :: iim_input,jjm_input,ngrid,nslope            ! number of points in the lat x lon dynamical grid, number of subgrid slopes
     27  INTEGER :: iim_input,jjm_input,ngrid,nslope      ! number of points in the lat x lon dynamical grid, number of subgrid slopes
    3128! Ouputs
    3229  REAL, INTENT(OUT) ::  min_co2_ice(ngrid,nslope) ! Minimum of co2 ice  per slope of the year [kg/m^2]
    … …  
    105102     print *, "Downloading data for surface pressure done"
    106103     print *, "nslope=", nslope
     104
     105if(nslope.gt.1) then
     106
    107107     print *, "Downloading data for co2ice_slope ..."
    108 
    109 if(nslope.gt.1) then
    110108
    111109DO islope=1,nslope
    … …  
    178176    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
    179177#ifndef CPP_STD
    180     call get_var3("watercap", watercap_slope(:,:,1,:))
     178!    call get_var3("watercap", watercap_slope(:,:,1,:))
     179  watercap_slope(:,:,1,:)=1.
    181180#endif
    182181
    … …  
    198197    tsurf_ave_dyn(:,:,:)=SUM(tsurf_gcm_dyn(:,:,:,:),4)/timelen
    199198
     199#ifndef CPP_1D
    200200  DO islope = 1,nslope
    201201    DO t=1,timelen
    … …  
    203203    ENDDO
    204204  ENDDO
     205#endif
    205206
    206207  if(soil_pem) then
    … …  
    213214! By definition, a density is positive, we get rid of the negative values
    214215  DO i=1,iim+1
    215     DO j = 1, jjm+1
     216    DO j = 1, jjm_input+1
    216217       DO islope=1,nslope
    217218          if (min_co2_ice_dyn(i,j,islope).LT.0) then
    … …  
    226227
    227228  DO i=1,iim+1
    228     DO j = 1, jjm+1
     229    DO j = 1, jjm_input+1
    229230      DO t = 1, timelen
    230231         if (q_co2_dyn(i,j,t).LT.0) then
    … …  
    243244    ENDDO
    244245  ENDDO
     246
     247#ifndef CPP_1D
    245248
    246249     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
    … …  
    269272     CALL gr_dyn_fi(nslope,iim_input+1,jjm_input+1,ngrid,tsurf_ave_dyn,tsurf_ave)
    270273
     274#else
     275
     276  vmr_co2_gcm_phys(1,:)=vmr_co2_gcm(1,1,:)
     277  ps_timeseries(1,:)=ps_GCM(1,1,:)
     278  q_co2(1,:)=q_co2_dyn(1,1,:)
     279  q_h2o(1,:)=q_h2o_dyn(1,1,:)
     280  min_co2_ice(1,:)=min_co2_ice_dyn(1,1,:)
     281  min_h2o_ice(1,:)=min_h2o_ice_dyn(1,1,:)
     282  if(soil_pem) then
     283    tsoil_ave(1,:,:)=tsoil_ave_dyn(1,1,:,:)
     284    tsoil_gcm(1,:,:,:)=tsoil_gcm_dyn(1,1,:,:,:)
     285    watersoil_density(1,:,:,:)=watersoil_density_dyn(1,1,:,:,:)
     286  endif !soil_pem
     287  tsurf_GCM(1,:,:)=tsurf_GCM_dyn(1,1,:,:)
     288  co2_ice_slope(1,:,:)=co2_ice_slope_dyn(1,1,:,:)
     289  tsurf_ave(1,:)=tsurf_ave_dyn(1,1,:)
     290
     291#endif
     292
    271293  CONTAINS
    272294

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

    r2963 r2980  
    33!
    44! Purpose: Recompute orbit parameters based on values by Laskar et al.,
    5 ! 2004
    65!
    7 ! Authors: RV
     6! Author: RV
    87!=======================================================================
    98      IMPLICIT NONE
    … …  
    4847      real :: unitastr
    4948
     49      real xa,xb,xc,ya,yb
     50
    5051
    5152      ! **********************************************************************
    … …  
    6667          print *, "recomp_orb_param, Old timeperi=", timeperi
    6768
     69        xc=Year
     70
    6871        do ilask=last_ilask,1,-1
     72
     73          xa=yearlask(ilask)
     74          xb=yearlask(ilask+1)
     75
    6976           if(yearlask(ilask) .GT.Year) then
    7077             if(var_obl) then
    71                obliquit=oblask(ilask+1)+(oblask(ilask)-oblask(ilask+1))*(Year-yearlask(ilask+1))/(yearlask(ilask)-yearlask(ilask+1))
     78              ya=oblask(ilask)
     79              yb=oblask(ilask+1)
     80              obliquit=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb)
    7281             endif
    7382             if(var_ex) then
    74                e_elips=exlask(ilask+1)+(exlask(ilask)-exlask(ilask+1))*(Year-yearlask(ilask+1))/(yearlask(ilask)-yearlask(ilask+1))
     83              ya=exlask(ilask)
     84              yb=exlask(ilask+1)
     85              e_elips=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb)
    7586             endif
    7687             if(var_lsp) then
    … …  
    8293                 endif
    8394               endif
    84                timeperi_ls=lsplask(ilask+1)+(lsplask(ilask)-lsplask(ilask+1))*(Year-yearlask(ilask+1))/(yearlask(ilask)-yearlask(ilask+1))
    85                timeperi_ls=modulo(timeperi_ls,360.)
     95
     96              ya=lsplask(ilask)
     97              yb=lsplask(ilask+1)
     98              timeperi_ls=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb)
     99              timeperi_ls=modulo(timeperi_ls,360.)
     100              halfaxe=227.94
     101              timeperi=timeperi_ls*2*pi/360
     102              periheli = halfaxe*(1-e_elips)
     103              aphelie =  halfaxe*(1+e_elips)
     104              unitastr=149.597927
     105              p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
     106
     107              call call_dayperi(timeperi,e_elips,peri_day,year_day)
    86108             endif
    87109              exit
    88110           endif
    89111        enddo
    90         halfaxe=227.94
    91         timeperi=timeperi_ls*2*pi/360
    92         periheli = halfaxe*(1-e_elips)
    93         aphelie =  halfaxe*(1+e_elips)
    94         unitastr=149.597927
    95         p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
    96         call call_dayperi(timeperi,e_elips,peri_day,year_day)
    97112
    98113       print *, "recomp_orb_param, Final year of the PEM run:", Year

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

    r2945 r2980  
    2727  REAL, intent(in) ::  tendencies_co2_ice_phys_ini(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
    2828  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)
     29  REAL, intent(in) :: emissivity_slope(ngrid,nslope)             ! Emissivity per mesh and sub-grid slope(1)
    3030
    3131!   OUTPUT
    … …  
    4343      coef=sols_per_my*sec_per_sol*emissivity_slope(i,islope)*sigmaB/Lco2
    4444      ave=0.
    45 !      if(abs(tendencies_co2_ice_phys(i,islope)).gt.1e-4) then
    4645      if(co2ice_slope(i,islope).gt.1e-4 .and. abs(tendencies_co2_ice_phys(i,islope)).gt.1e-5) then
    4746        do t=1,timelen

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

    r2963 r2980  
    8080  status = nf90_inquire_dimension(ncid1, dimids(i), name_, len_)
    8181  if(status /= nf90_noerr) call handle_err(status)
    82   if(name_.eq."lon")  then
     82  if(name_.eq."lon" .or. name_.eq."longitude")  then
    8383     dimid_lon=dimids(i)
    8484     len_lon=len_
    8585     len_=len_+1
    86   elseif(name_.eq."lat") then
     86  elseif(name_.eq."lat".or. name_.eq."latitude") then
    8787     dimid_lat=dimids(i)
    8888     len_lat=len_
    89   elseif(name_.eq."time_counter") then
     89  elseif(name_.eq."time_counter".or. name_.eq. "Time") then
    9090     dimid_time=dimids(i)
    9191     len_time=len_
    92   elseif(name_.eq."soil_layers") then
     92  elseif(name_.eq."soil_layers".or. name_.eq. "subsurface_layers") then
    9393     dimid_soil=dimids(i)
    9494     len_soil=len_
    … …  
    9898  dimids_2(i)=dimid_2
    9999enddo
    100 
    101 
    102 
    103 allocate(tempvalues_3d(len_lon,len_lat,len_time))
    104 allocate(values_3d(len_lon+1,len_lat,len_time))
    105 
    106 allocate(tempvalues_4d(len_lon,len_lat,len_soil,len_time))
    107 allocate(values_4d(len_lon+1,len_lat,len_soil,len_time))
    108 
    109100
    110101do i=1,nvars
    … …  
    186177    endif
    187178  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))
    188181      status = nf90_get_var(ncid1, varids(i), tempvalues_3d)
    189182      if(status /= nf90_noerr) call handle_err(status)
    … …  
    198191      status = nf90_redef(ncid2)
    199192      if(status /= nf90_noerr) call handle_err(status)
     193      deallocate(tempvalues_3d)
     194      deallocate(values_3d)
    200195  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))
    201198      status = nf90_get_var(ncid1, varids(i), tempvalues_4d)
    202199      if(status /= nf90_noerr) call handle_err(status)
    … …  
    211208      status = nf90_redef(ncid2)
    212209      if(status /= nf90_noerr) call handle_err(status)
     210      deallocate(tempvalues_4d)
     211      deallocate(values_4d)
    213212  endif
    214213
    … …  
    228227deallocate(dimids_2)
    229228deallocate(varids_2)
    230 deallocate(tempvalues_3d)
    231 deallocate(values_3d)
    232 deallocate(tempvalues_4d)
    233 deallocate(values_4d)
    234 
    235 
    236229
    237230enddo
Note: See TracChangeset for help on using the changeset viewer.