Changeset 3591


Ignore:
Timestamp:
Jan 21, 2025, 3:53:18 PM (15 hours ago)
Author:
jbclement
Message:

PEM:

  • Making allocation/deallocation systematically and more efficient in the main program.
  • Some cleanings (variables deletion, more adapted type/dimension, etc).

JBC

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

Legend:

Unmodified
Added
Removed

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

    r3589 r3591  
    547547== 20/01/2025 == JBC
    548548The initialization of constants between dynamics and physics was messed up which led to uninitialized variables.
     549
     550== 21/01/2025 == JBC
     551- Making allocation/deallocation systematically and more efficient in the main program.
     552- Some cleanings (variables deletion, more adapted type/dimension, etc).

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

    r3571 r3591  
    4848        write(*,*)"read_data_PCM: Le champ <time_counter> est absent"
    4949        write(*,*)"read_data_PCM: J essaie <Time>"
    50         ierr = nf90_inq_varid (fID, "Time", vID)
     50        ierr = nf90_inq_varid (fID,"Time",vID)
    5151        if (ierr /= nf90_noerr) then
    5252            write(*,*)"read_data_PCM: Le champ <Time> est absent"

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

    r3571 r3591  
    2121!=======================================================================
    2222
    23 SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,co2ice,flag_co2flow,flag_co2flow_mesh)
     23SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,co2ice,flag_co2flow)
    2424
    2525!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    3535implicit none
    3636
    37 ! Inputs:
    38 !--------
     37! Inputs
     38!-------
    3939integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
    4040real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Grid points x Slopes: Distribution of the subgrid slopes
    … …  
    4444real,                              intent(in) :: ps_avg_global_ini             ! Global averaged surface pressure at the beginning [Pa]
    4545real,                              intent(in) :: ps_avg_global                 ! Global averaged surface pressure during the PEM iteration [Pa]
    46 
    47 ! Ouputs:
    48 !--------
    49 real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    50 real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
    51 real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
     46! Ouputs
     47!-------
     48real, dimension(ngrid,nslope), intent(inout) :: co2ice ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     49integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_co2flow ! Flag to see if there is flow on the subgrid slopes
    5250! Local
    5351!------
    … …  
    5856call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
    5957call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
    60 call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tcond,co2ice,flag_co2flow,flag_co2flow_mesh)
     58call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tcond,co2ice,flag_co2flow)
    6159
    6260END SUBROUTINE flow_co2glaciers
    … …  
    6462!=======================================================================
    6563
    66 SUBROUTINE flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
     64SUBROUTINE flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow)
    6765
    6866!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    8179! ---------
    8280
    83 ! Inputs:
     81! Inputs
     82! ------
    8483integer,                       intent(in) :: ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
    8584real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Grid points x Slopes : Distribution of the subgrid slopes
    8685real, dimension(ngrid),        intent(in) :: def_slope_mean ! Slopes: values of the sub grid slope angles
    8786real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice Temperature [K]
    88 ! Ouputs:
    89 real, dimension(ngrid,nslope), intent(inout) :: h2oice            ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    90 real, dimension(ngrid,nslope), intent(inout) :: flag_h2oflow      ! flag to see if there is flow on the subgrid slopes
    91 real, dimension(ngrid),        intent(inout) :: flag_h2oflow_mesh ! same but within the mesh
    92 ! Local
     87! Outputs
     88!--------
     89real, dimension(ngrid,nslope), intent(inout) :: h2oice ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     90integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_h2oflow ! Flag to see if there is flow on the subgrid slopes
     91! Local
     92! -----
    9393real, dimension(ngrid,nslope) :: hmax ! Grid points x Slope field: maximum thickness for co2  glacier before flow
    9494
    9595write(*,*) "Flow of H2O glaciers"
    9696call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
    97 call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
     97call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,h2oice,flag_h2oflow)
    9898
    9999END SUBROUTINE flow_h2oglaciers
    … …  
    115115!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle before initating flow
    116116!!!
    117 !!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
     117!!! Author: LL, based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
    118118!!!
    119119!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    161161!=======================================================================
    162162
    163 SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,qice,flag_flow,flag_flowmesh)
     163SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,qice,flag_flow)
    164164!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    165165!!!
    … …  
    190190! Outputs
    191191!--------
    192 real, dimension(ngrid,nslope), intent(inout) :: qice          ! CO2 in the subslope [kg/m^2]
    193 real, dimension(ngrid,nslope), intent(inout) :: flag_flow     ! boolean to check if there is flow on a subgrid slope
    194 real, dimension(ngrid),        intent(inout) :: flag_flowmesh ! boolean to check if there is flow in the mesh
     192real, dimension(ngrid,nslope), intent(inout) :: qice ! CO2 in the subslope [kg/m^2]
     193integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_flow ! Flag to see if there is flow on the subgrid slopes
    195194! Local
    196195!------
    197196integer :: ig, islope ! loop
    198197integer :: iaval      ! ice will be transfered here
     198
     199flag_flow = 0
    199200
    200201do ig = 1,ngrid
    … …  
    220221
    221222                qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)
    222                 flag_flow(ig,islope) = 1.
    223                 flag_flowmesh(ig) = 1.
     223                flag_flow(ig,islope) = 1
    224224            endif ! co2ice > hmax
    225225        endif ! iflat

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

    r3589 r3591  
    179179real,   dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Grid points x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
    180180integer                               :: timelen              ! # time samples
    181 real,   dimension(:,:),   allocatable :: p                    ! Grid points x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
    182181real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
    183182
    … …  
    200199
    201200! Variables for slopes
    202 real, dimension(:,:),   allocatable :: flag_co2flow      ! (ngrid,nslope): Flag where there is a CO2 glacier flow
    203 real, dimension(:),     allocatable :: flag_co2flow_mesh ! (ngrid)       : Flag where there is a CO2 glacier flow
    204 real, dimension(:,:),   allocatable :: flag_h2oflow      ! (ngrid,nslope): Flag where there is a H2O glacier flow
    205 real, dimension(:),     allocatable :: flag_h2oflow_mesh ! (ngrid)       : Flag where there is a H2O glacier flow
     201integer(kind = 1), dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow
     202integer(kind = 1), dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow
    206203
    207204! Variables for surface and soil
    … …  
    218215real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
    219216real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Grid points x Soil x Slopes: Averaged water soil density [kg/m^3]
    220 real, dimension(:,:),     allocatable :: tsurf_avg_old                    ! Surface temperature saved from previous time step [K]
    221217real, dimension(:),       allocatable :: delta_co2_adsorbed               ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    222218real, dimension(:),       allocatable :: delta_h2o_adsorbed               ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
    … …  
    233229
    234230! Some variables for the PEM run
    235 real, parameter :: year_step = 1   ! Timestep for the pem
    236 real            :: i_myear_leg     ! Number of iteration
    237 real            :: n_myear_leg     ! Maximum number of iterations before stopping
    238 real            :: i_myear         ! Global number of Martian years of the chained simulations
    239 real            :: n_myear         ! Maximum number of Martian years of the chained simulations
    240 real            :: timestep        ! Timestep [s]
    241 character(100)  :: arg             ! To read command-line arguments program was invoked
    242 logical         :: timewall        ! Flag to use the time limit stopping criterion in case of a PEM job
    243 integer(kind=8) :: cr              ! Number of clock ticks per second (count rate)
    244 integer(kind=8) :: c1, c2          ! Counts of processor clock
    245 character(100)  :: chtimelimit     ! Time limit for the PEM job outputted by the SLURM command
    246 real            :: timelimit       ! Time limit for the PEM job in seconds
    247 real, parameter :: antetime = 1200 ! Anticipation time to prevent reaching the time limit: 1200 s = 20 min by default
    248 integer         :: cstat, days, hours, minutes, seconds
    249 character(1)    :: sep
     231real, parameter   :: year_step = 1   ! Timestep for the pem
     232real              :: i_myear_leg     ! Number of iteration
     233real              :: n_myear_leg     ! Maximum number of iterations before stopping
     234real              :: i_myear         ! Global number of Martian years of the chained simulations
     235real              :: n_myear         ! Maximum number of Martian years of the chained simulations
     236real              :: timestep        ! Timestep [s]
     237character(100)    :: arg             ! To read command-line arguments program was invoked
     238logical           :: timewall        ! Flag to use the time limit stopping criterion in case of a PEM job
     239integer(kind = 8) :: cr              ! Number of clock ticks per second (count rate)
     240integer(kind = 8) :: c1, c2          ! Counts of processor clock
     241character(100)    :: chtimelimit     ! Time limit for the PEM job outputted by the SLURM command
     242real              :: timelimit       ! Time limit for the PEM job in seconds
     243real, parameter   :: antetime = 1200 ! Anticipation time to prevent reaching the time limit: 1200 s = 20 min by default
     244integer           :: cstat, days, hours, minutes, seconds
     245character(1)      :: sep
    250246
    251247#ifdef CPP_STD
    … …  
    281277    real, dimension(nlayer + 1)         :: plev
    282278#else
    283     integer, parameter              :: jjm_value = jjm
    284     real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart
    285     real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart
     279    integer, parameter                :: jjm_value = jjm
     280    real, dimension(:),   allocatable :: ap ! Coefficient ap read in start_name and written in restart
     281    real, dimension(:),   allocatable :: bp ! Coefficient bp read in start_name and written in restart
     282    real, dimension(:,:), allocatable :: p  ! Grid points x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
    286283#endif
    287284
    … …  
    346343
    347344! Some constants
    348 day_ini = 0    ! test
    349 time_phys = 0. ! test
     345day_ini = 0
     346time_phys = 0.
    350347ngrid = ngridmx
    351 A = (1/m_co2 - 1/m_noco2)
    352 B = 1/m_noco2
     348A = (1./m_co2 - 1./m_noco2)
     349B = 1./m_noco2
    353350year_day = 669
    354351daysec = 88775.
    … …  
    399396!------------------------
    400397! I_b.1 Read "start.nc"
    401 allocate(ps_start(ngrid),ps_start0(ngrid))
     398allocate(ps_start0(ngrid))
    402399#ifndef CPP_1D
    403400    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps_start_dyn,phis,time_0)
    404401
    405     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps_start_dyn,ps_start)
     402    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps_start_dyn,ps_start0)
    406403
    407404    call iniconst ! Initialization of dynamical constants (comconst_mod)
    … …  
    417414    status = nf90_close(ncid)
    418415
    419     ! Initialization of physics constants (comcstfi_h)
     416    ! Initialization of physics constants and variables (comcstfi_h)
    420417    call iniphysiq(iim,jjm,llm,(jjm - 1)*iim + 2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
    421418#else
    422     ps_start(1) = ps_start_dyn(1)
     419    ps_start0(1) = ps_start_dyn(1)
    423420#endif
    424 
    425 ! Save initial surface pressure
    426 ps_start0 = ps_start
    427421
    428422! In the PCM, these values are given to the physic by the dynamic.
    429423! Here we simply read them in the "startfi.nc" file
    430424status = nf90_open(startfi_name,NF90_NOWRITE,ncid)
    431 
    432425status = nf90_inq_varid(ncid,"longitude",lonvarid)
    433426status = nf90_get_var(ncid,lonvarid,longitude)
    434 
    435427status = nf90_inq_varid(ncid,"latitude",latvarid)
    436428status = nf90_get_var(ncid,latvarid,latitude)
    437 
    438429status = nf90_inq_varid(ncid,"area",areavarid)
    439430status = nf90_get_var(ncid,areavarid,cell_area)
    440 
    441431status = nf90_inq_varid(ncid,"soildepth",sdvarid)
    442432status = nf90_get_var(ncid,sdvarid,mlayer)
    443 
    444433status = nf90_close(ncid)
    445434
    … …  
    465454    allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
    466455    allocate(emis_read_generic(ngrid))
     456    allocate(albedo_read_generic(ngrid,2))
     457    allocate(qsurf(ngrid,nq,1))
    467458    allocate(tsurf(ngrid,1))
    468     allocate(qsurf(ngrid,nq,1))
    469459    allocate(tsoil(ngrid,nsoilmx,1))
    470460    allocate(emis(ngrid,1))
    471461    allocate(watercap(ngrid,1))
    472462    allocate(watercaptag(ngrid))
    473     allocate(albedo_read_generic(ngrid,2))
    474463    allocate(albedo(ngrid,2,1))
    475464    allocate(inertiesoil(ngrid,nsoilmx,1))
    … …  
    503492    qsurf(:,:,1) = qsurf_read_generic
    504493    where (qsurf < 0.) qsurf = 0.
     494
     495    deallocate(tsurf_read_generic,qsurf_read_generic,qsoil_read_generic,emis_read_generic)
    505496#endif
    506497
    … …  
    523514    if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
    524515enddo
    525 
    526516write(*,*) 'Flat slope for islope = ',iflat
    527517write(*,*) 'corresponding criterium = ',def_slope_mean(iflat)
    528518
    529 allocate(flag_co2flow(ngrid,nslope))
    530 allocate(flag_co2flow_mesh(ngrid))
    531 allocate(flag_h2oflow(ngrid,nslope))
    532 allocate(flag_h2oflow_mesh(ngrid))
    533 
    534 flag_co2flow = 0
    535 flag_co2flow_mesh = 0
    536 flag_h2oflow = 0
    537 flag_h2oflow_mesh = 0
    538 
    539519!------------------------
    540520! I   Initialization
    … …  
    544524call get_timelen_PCM("data_PCM_Y1.nc",timelen)
    545525
    546 allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope))
    547 allocate(vmr_co2_PCM(ngrid,timelen))
    548 allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid),ps_dev(ngrid))
     526allocate(vmr_co2_PCM(ngrid,timelen),q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
     527allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid))
    549528allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
    550 allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope),tsurf_avg_old(ngrid,nslope),tsurf_dev(ngrid,nslope))
    551 allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_dev(ngrid,nsoilmx,nslope))
    552 allocate(tsoil_timeseries(ngrid,nsoilmx,nslope,timelen),tsoil_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    553 allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
    554 allocate(q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
    555 allocate(watersurf_density_avg(ngrid,nslope))
    556 allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen),watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    557 allocate(delta_co2_adsorbed(ngrid))
    558 allocate(co2ice_disappeared(ngrid,nslope))
    559 allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
    560 allocate(delta_h2o_icetablesublim(ngrid),delta_h2o_adsorbed(ngrid))
    561 allocate(vmr_co2_PEM_phys(ngrid,timelen))
     529allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope))
     530allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_timeseries(ngrid,nsoilmx,nslope,timelen),tsoil_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
     531allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
    562532
    563533call read_data_PCM("data_PCM_Y1.nc","data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,ps_avg,tsurf_avg_yr1,tsurf_avg, &
    … …  
    565535
    566536! Compute the deviation from the average
    567 ps_dev = ps_start - ps_avg
     537allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoilmx,nslope))
     538ps_dev = ps_start0 - ps_avg
    568539tsurf_dev = tsurf - tsurf_avg
    569540tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:)
    … …  
    580551call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
    581552
     553allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope),watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    582554if (soil_pem) then
    583555    call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
    … …  
    598570!    I_f Compute tendencies
    599571!------------------------
    600 allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope))
    601 allocate(d_co2ice_ini(ngrid,nslope))
    602 
    603 ! Compute the tendencies of the evolution of ice over the years
     572allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope))
    604573call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
    605574call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
    … …  
    621590!    I_h Read the "startpem.nc"
    622591!------------------------
    623 allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
    624 allocate(stratif(ngrid,nslope))
     592allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),stratif(ngrid,nslope))
     593allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid))
    625594if (layering_algo) then
    626595    do islope = 1,nslope
    … …  
    631600endif
    632601
     602delta_h2o_icetablesublim = 0.
    633603call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
    634604              icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
    … …  
    639609! We save the places where h2o ice is sublimating
    640610! We compute the surface of h2o ice sublimating
    641 allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
     611allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope),co2ice_disappeared(ngrid,nslope))
    642612co2ice_sublim_surf_ini = 0.
    643613h2oice_ini_surf = 0.
    … …  
    662632write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
    663633
    664 delta_h2o_icetablesublim = 0.
    665634
    666635if (adsorption_pem) then
    … …  
    763732! II.a.3. Tracers timeseries
    764733    write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..."
     734    allocate(vmr_co2_PEM_phys(ngrid,timelen))
    765735    l = 1
    766736    do ig = 1,ngrid
    … …  
    795765!    II_b Evolution of ice
    796766!------------------------
     767    allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
    797768    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
    798769    call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
    … …  
    811782!    II_c Flow of glaciers
    812783!------------------------
     784    allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope))
    813785    if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
    814                                                             ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
    815     if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
     786                                                            ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
     787    if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
    816788
    817789!------------------------
    … …  
    854826! II_d.2 Shifting soil temperature to surface
    855827        call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
     828        deallocate(zshift_surf,zlag)
    856829
    857830! II_d.3 Update soil temperature
    … …  
    880853
    881854! II_d.4 Update the ice table
     855        allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
    882856        if (icetable_equilibrium) then
    883857            write(*,*) "Compute ice table (equilibrium method)"
    … …  
    899873            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
    900874        endif
     875        deallocate(icetable_thickness_old,ice_porefilling_old)
    901876
    902877! II_d.5 Update the soil thermal properties
    … …  
    944919        call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
    945920        call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
    946         call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
     921        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope)))
     922        call writediagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope)))
    947923        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
    948924        if (icetable_equilibrium) then
    … …  
    963939        endif
    964940    enddo
     941    deallocate(flag_co2flow,flag_h2oflow)
    965942
    966943!------------------------
    … …  
    969946!------------------------
    970947    call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_old,ps_avg_global_new)
     948    deallocate(vmr_co2_PEM_phys)
    971949
    972950!------------------------
    … …  
    1011989    endif
    1012990enddo ! big time iteration loop of the pem
     991deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed)
     992deallocate(watersoil_density_PEM_avg,watersurf_density_avg,)
     993deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries)
     994deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim)
     995deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
     996deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
     997if (layering_algo) then
     998    do islope = 1,nslope
     999        do i = 1,ngrid
     1000            call del_layering(stratif(i,islope))
     1001        enddo
     1002    enddo
     1003    deallocate(new_str,new_lag1,new_lag2,current1,current2)
     1004endif
    10131005!------------------------------ END RUN --------------------------------
    10141006
    … …  
    10481040! III.a.2. Tsurf update for start file
    10491041tsurf = tsurf_avg + tsurf_dev
     1042deallocate(tsurf_dev)
    10501043
    10511044! III_a.3 Tsoil update for start file
    … …  
    10611054#endif
    10621055endif
     1056deallocate(tsurf_avg,tsoil_dev)
    10631057
    10641058! III_a.4 Pressure update for start file
     1059allocate(ps_start(ngrid))
    10651060ps_start = ps_avg + ps_dev
     1061deallocate(ps_avg,ps_dev)
    10661062
    10671063! III_a.5 Tracers update for start file
    … …  
    10901086    endif
    10911087enddo
     1088deallocate(zplev_start0)
    10921089
    10931090! Conserving the tracers mass for start file
    … …  
    11051102    enddo
    11061103enddo
    1107 deallocate(zplev_start0,zplev_new)
     1104deallocate(zplev_new)
    11081105
    11091106! III_a.6 Albedo update for start file
    … …  
    11551152ztime_fin = time_phys
    11561153
    1157 allocate(p(ip1jmp1,nlayer + 1))
    11581154#ifndef CPP_1D
     1155    allocate(p(ip1jmp1,nlayer + 1))
    11591156    ! Correction on teta due to surface pressure changes
    11601157    do l = 1,nlayer
    … …  
    11691166    call dynredem0("restart.nc",day_ini,phis)
    11701167    call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps_start)
    1171     write(*,*) "restart.nc has been written"
     1168    write(*,*) "restart.nc has been written."
     1169    deallocate(ap,bp,p)
    11721170#else
    11731171    call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
    1174     write(*,*) "restart1D.txt has been written"
     1172    write(*,*) "restart1D.txt has been written."
    11751173#endif
     1174deallocate(ps_start0,ps_start)
    11761175
    11771176! III_b.2 Write the "restartfi.nc"
    … …  
    11851184                  wstar,watercap,perennial_co2ice)
    11861185#else
     1186    if (allocated(noms)) deallocate(noms)
     1187    deallocate(qsurf,tsurf,tsoil,emis,watercap,watercaptag,albedo,inertiesoil)
    11871188    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
    11881189                  nlayer,nq,ptimestep,pday,time_phys,cell_area,    &
    … …  
    11931194                  tsea_ice,sea_ice)
    11941195#endif
    1195 write(*,*) "restartfi.nc has been written"
     1196write(*,*) "restartfi.nc has been written."
    11961197
    11971198!------------------------
    … …  
    12031204call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
    12041205             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,stratif)
    1205 write(*,*) "restartpem.nc has been written"
     1206write(*,*) "restartpem.nc has been written."
    12061207
    12071208call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
    … …  
    12181219        enddo
    12191220    enddo
    1220     deallocate(new_str,new_lag1,new_lag2,current1,current2)
    12211221endif
    1222 deallocate(ps_start,ps_start0,ps_timeseries,ps_avg,ps_dev)
    1223 deallocate(tsurf_avg,tsurf_dev,tsurf_avg_old)
    1224 deallocate(tsoil_PEM,tsoil_dev,tsoil_PEM_timeseries)
    1225 deallocate(vmr_co2_PCM,vmr_co2_PEM_phys,q_co2_PEM_phys,q_h2o_PEM_phys)
    1226 deallocate(watersurf_density_avg,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
    1227 deallocate(delta_co2_adsorbed,delta_h2o_adsorbed,delta_h2o_icetablesublim)
    1228 deallocate(icetable_thickness_old,ice_porefilling_old,zshift_surf,zlag)
    1229 deallocate(is_co2ice_ini,co2ice_disappeared,is_co2ice_sublim_ini,is_h2oice_sublim_ini,stratif)
     1222deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
     1223deallocate(co2_ice,h2o_ice,stratif)
    12301224!----------------------------- END OUTPUT ------------------------------
    12311225

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

    r3571 r3591  
    121121        call put_field("mco2_reg_ads_slope"//num, "Mass of CO2 adsorbed in the regolith",m_co2_regolith(:,:,islope),Year)
    122122        call put_field("mh2o_reg_ads_slope"//num, "Mass of H2O adsorbed in the regolith",m_h2o_regolith(:,:,islope),Year)
     123        call put_field("ice_porefilling"//num,"Subsurface ice pore filling",ice_porefilling(:,:,islope),Year)
    123124    enddo
    124125    call put_field("icetable_depth","Depth of ice table",icetable_depth,Year)
    125126    call put_field("icetable_thickness","Depth of ice table",icetable_thickness,Year)
    126     call put_field("ice_porefilling","Subsurface ice pore filling",ice_porefilling,Year)
    127127    call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
    128128endif ! soil_pem

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

    r3571 r3591  
    8484if (nslope == 1) then ! There is no slope
    8585    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
    86     write(*,*) "Data for co2_ice downloaded!"
     86    write(*,*) "Data for co2_ice downloaded."
    8787
    8888    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
    89     write(*,*) "Data for h2o_ice_s downloaded!"
     89    write(*,*) "Data for h2o_ice_s downloaded."
    9090
    9191    call get_var3("tsurf",tsurf_dyn(:,:,1,:))
    92     write(*,*) "Data for tsurf downloaded!"
     92    write(*,*) "Data for tsurf downloaded."
    9393
    9494#ifndef CPP_STD
    9595    call get_var3("watercap",watercap(:,:,1,:))
    96     write(*,*) "Data for watercap downloaded!"
     96    write(*,*) "Data for watercap downloaded."
    9797
    9898    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
    99     write(*,*) "Data for perennial_co2ice downloaded!"
     99    write(*,*) "Data for perennial_co2ice downloaded."
    100100#endif
    101101else ! We use slopes
    … …  
    104104        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
    105105    enddo
    106     write(*,*) "Data for co2_ice downloaded!"
     106    write(*,*) "Data for co2_ice downloaded."
    107107
    108108    do islope = 1,nslope
    … …  
    110110        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
    111111    enddo
    112     write(*,*) "Data for h2o_ice_s downloaded!"
     112    write(*,*) "Data for h2o_ice_s downloaded."
    113113
    114114    do islope = 1,nslope
    … …  
    116116        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
    117117    enddo
    118     write(*,*) "Data for tsurf downloaded!"
     118    write(*,*) "Data for tsurf downloaded."
    119119
    120120#ifndef CPP_STD
    … …  
    123123        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
    124124    enddo
    125     write(*,*) "Data for watercap downloaded!"
     125    write(*,*) "Data for watercap downloaded."
    126126
    127127    do islope = 1,nslope
    … …  
    129129        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
    130130    enddo
    131     write(*,*) "Data for perennial_co2ice downloaded!"
     131    write(*,*) "Data for perennial_co2ice downloaded."
    132132#endif
    133133endif
    … …  
    170170! Download the data from the file
    171171call get_var3("ps",ps_dyn)
    172 write(*,*) "Data for surface pressure downloaded!"
     172write(*,*) "Data for surface pressure downloaded."
    173173
    174174call get_var3("co2_layer1",q_co2_dyn)
    175 write(*,*) "Data for vmr co2 downloaded!"
     175write(*,*) "Data for vmr co2 downloaded."
    176176
    177177call get_var3("h2o_layer1",q_h2o_dyn)
    178 write(*,*) "Data for vmr h2o downloaded!"
     178write(*,*) "Data for vmr h2o downloaded."
    179179
    180180if (nslope == 1) then ! There is no slope
    181181    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
    182     write(*,*) "Data for co2_ice downloaded!"
     182    write(*,*) "Data for co2_ice downloaded."
    183183
    184184    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
    185     write(*,*) "Data for h2o_ice_s downloaded!"
     185    write(*,*) "Data for h2o_ice_s downloaded."
    186186
    187187    call get_var3("tsurf",tsurf_dyn(:,:,1,:))
    188     write(*,*) "Data for tsurf downloaded!"
     188    write(*,*) "Data for tsurf downloaded."
    189189
    190190#ifndef CPP_STD
    191191    call get_var3("watercap",watercap(:,:,1,:))
    192     write(*,*) "Data for watercap downloaded!"
     192    write(*,*) "Data for watercap downloaded."
    193193
    194194    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
    195     write(*,*) "Data for perennial_co2ice downloaded!"
     195    write(*,*) "Data for perennial_co2ice downloaded."
    196196
    197197    if (soil_pem) then
    198198        call get_var4("soiltemp",tsoil_dyn(:,:,:,1,:))
    199         write(*,*) "Data for soiltemp downloaded!"
     199        write(*,*) "Data for soiltemp downloaded."
    200200
    201201        call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
    202         write(*,*) "Data for waterdensity_soil downloaded!"
     202        write(*,*) "Data for waterdensity_soil downloaded."
    203203
    204204        call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,1,:))
    205         write(*,*) "Data for waterdensity_surface downloaded!"
     205        write(*,*) "Data for waterdensity_surface downloaded."
    206206    endif !soil_pem
    207207#endif
    … …  
    211211        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
    212212    enddo
    213     write(*,*) "Data for co2_ice downloaded!"
     213    write(*,*) "Data for co2_ice downloaded."
    214214
    215215    do islope = 1,nslope
    … …  
    217217        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
    218218    enddo
    219     write(*,*) "Data for h2o_ice_s downloaded!"
     219    write(*,*) "Data for h2o_ice_s downloaded."
    220220
    221221    do islope = 1,nslope
    … …  
    223223        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
    224224    enddo
    225     write(*,*) "Data for tsurf downloaded!"
     225    write(*,*) "Data for tsurf downloaded."
    226226
    227227#ifndef CPP_STD
    … …  
    230230        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
    231231    enddo
    232     write(*,*) "Data for watercap downloaded!"
     232    write(*,*) "Data for watercap downloaded."
    233233
    234234    do islope = 1,nslope
    … …  
    236236        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
    237237    enddo
    238     write(*,*) "Data for perennial_co2ice downloaded!"
     238    write(*,*) "Data for perennial_co2ice downloaded."
    239239
    240240    if (soil_pem) then
    … …  
    243243            call get_var4("soiltemp_slope"//num,tsoil_dyn(:,:,:,islope,:))
    244244        enddo
    245         write(*,*) "Data for soiltemp downloaded!"
     245        write(*,*) "Data for soiltemp downloaded."
    246246
    247247        do islope = 1,nslope
    … …  
    249249            call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
    250250        enddo
    251         write(*,*) "Data for waterdensity_soil downloaded!"
     251        write(*,*) "Data for waterdensity_soil downloaded."
    252252
    253253        do islope = 1,nslope
    … …  
    255255            call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
    256256        enddo
    257         write(*,*) "Data for waterdensity_surface downloaded!"
     257        write(*,*) "Data for waterdensity_surface downloaded."
    258258    endif !soil_pem
    259259#endif
Note: See TracChangeset for help on using the changeset viewer.