Changeset 2963


Ignore:
Timestamp:
May 11, 2023, 5:40:03 PM (20 months ago)
Author:
romain.vande
Message:

Mars PEM :

Adapt PEM to the subslope PCM configuration, it is now fully compatible.

Create a PEM folder in deftank that contains:

run_pem1: a bash file that runs chained simulation of PEM as well as running a parameterizable number of PCM simulation in between.

It also takes care of reshaping XIOS output as well as renaming outputs etc… in the spirit of run_month1.

It is written for Irene machine and the header needs to be adapted for other machines.

run_PEM.def: A text file that shows the possible parameters to choose before a PEM simulation.

It should be included at the end of run.def just like callphys.def

ob_ex_lsp.asc: An ascii file containing the obliquity, eccentricity, ls_peri data from Laskar in Martian year.
README: A txt file explaining the content of the folder

Adapt field_def_physics_mars.xml to consider the case with 7 subslopes in the PCM.
Change context_lmdz_physics.xml to be able to output the file needed by the PEM.

Correct a few other minor bugs.

RV

Location:
trunk
Files:
5 added
14 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90

    r2961 r2963  
    2525!PEM parameters
    2626
    27 
    28 
    29   year_bp_ini=0.
    30   CALL getin('year_bp_ini', year_bp_ini)
    31 
    32   dt_pem=1
    33   CALL getin('dt_pem', dt_pem)
    34 
    35   water_ice_criterion=0.2
    36   CALL getin('water_ice_criterion', water_ice_criterion)
    37 
    38   co2_ice_criterion=0.2
    39   CALL getin('co2_ice_criterion', co2_ice_criterion)
    40 
    41   ps_criterion = 0.15
    42   CALL getin('ps_criterion',ps_criterion)
     27! #---------- ORBITAL parameters --------------#
    4328
    4429  evol_orbit_pem=.false.
    4530  CALL getin('evol_orbit_pem', evol_orbit_pem)
    4631
    47   Max_iter_pem=99999999
    48   CALL getin('Max_iter_pem', Max_iter_pem)
    49 
    50   co2glaciersflow = .true.
    51   CALL getin('co2glaciersflow', co2glaciersflow)
    52 
    53   soil_pem=.true.
    54   CALL getin('soil_pem', soil_pem)
    55 
    56   adsorption_pem = .true.
    57   CALL getin('adsorption_pem',adsorption_pem)
    58 
    59   fluxgeo = 0.
    60   CALL getin('Fluxgeo_PEM',fluxgeo)
    61   print*,'Flux Geothermal is set to',fluxgeo
     32  year_bp_ini=0.
     33  CALL getin('year_bp_ini', year_bp_ini)
    6234
    6335  var_obl = .true.
     
    7244  CALL getin('var_lsp',var_lsp)
    7345  print*,'Does Ls peri vary ?',var_lsp
     46
     47! #---------- Stopping criterion parameters --------------#
     48
     49  Max_iter_pem=99999999
     50  CALL getin('Max_iter_pem', Max_iter_pem)
     51
     52  water_ice_criterion=0.2
     53  CALL getin('water_ice_criterion', water_ice_criterion)
     54
     55  co2_ice_criterion=0.2
     56  CALL getin('co2_ice_criterion', co2_ice_criterion)
     57
     58  ps_criterion = 0.15
     59  CALL getin('ps_criterion',ps_criterion)
     60
     61  dt_pem=1
     62  CALL getin('dt_pem', dt_pem)
     63
     64! #---------- Subsurface parameters --------------#
     65
     66  soil_pem=.true.
     67  CALL getin('soil_pem', soil_pem)
     68
     69  adsorption_pem = .true.
     70  CALL getin('adsorption_pem',adsorption_pem)
     71
     72  co2glaciersflow = .true.
     73  CALL getin('co2glaciersflow', co2glaciersflow)
     74
     75  reg_thprop_dependp = .false.
     76  CALL getin('reg_thprop_dependp',reg_thprop_dependp)
     77  print*, 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
     78
     79! #---------- Layering parameters --------------#
     80
     81  fluxgeo = 0.
     82  CALL getin('Fluxgeo_PEM',fluxgeo)
     83  print*,'Flux Geothermal is set to',fluxgeo
    7484   
    7585  depth_breccia   = 10.
     
    8090  CALL getin('depth_bedrock',depth_bedrock)
    8191  print*,'Depth of bedrock is set to',depth_bedrock
    82 
    83   reg_thprop_dependp = .false.
    84   CALL getin('reg_thprop_dependp',reg_thprop_dependp)
    85   print*, 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
    8692
    8793   icetable_equilibrium = .true.
  • trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90

    r2950 r2963  
    8686    enddo
    8787   endif
    88   endif
    8988  negative_part = 0.
    9089
  • trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90

    r2894 r2963  
    4646
    4747      logical obl_not_found, ex_not_found,lsp_not_found !Loop variable (first call)
     48      logical max_lsp_modulo,min_lsp_modulo
    4849
    4950      ! **********************************************************************
     
    109110        min_lsp=timeperi_ls-max_change_lsp
    110111
     112        max_lsp_modulo=.false.
     113        min_lsp_modulo=.false.
     114        if(max_lsp.gt.360) max_lsp_modulo=.true.
     115        if(min_lsp.lt.0) min_lsp_modulo=.true.
     116
    111117!End Constant max change case
    112118
     
    133139        max_lsp_iter=999999999999
    134140
    135         do ilask=last_ilask+1,1,-1
     141        do ilask=last_ilask,1,-1
    136142           if((oblask(ilask).GT.max_obl).and. obl_not_found ) then
    137               max_obl_iter=(max_obl-oblask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
    138                                / (oblask(ilask-1)-oblask(ilask))
     143              max_obl_iter=(max_obl-oblask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
     144                               / (oblask(ilask)-oblask(ilask-1))+yearlask(ilask-1)-Year
    139145              obl_not_found=.FALSE.
    140146           elseif((oblask(ilask).LT.min_obl).and. obl_not_found ) then
    141               max_obl_iter=(min_obl-oblask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
    142                                / (oblask(ilask-1)-oblask(ilask))
     147              max_obl_iter=(min_obl-oblask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
     148                               / (oblask(ilask)-oblask(ilask-1))+yearlask(ilask-1)-Year
    143149              obl_not_found=.FALSE.
    144150           endif
    145151           if((exlask(ilask).GT.max_ex).and. ex_not_found ) then
    146               max_ex_iter=(max_ex-exlask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
    147                                / (exlask(ilask-1)-exlask(ilask))
     152              max_ex_iter=(max_ex-exlask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
     153                               / (exlask(ilask)-exlask(ilask-1))+yearlask(ilask-1)-Year
    148154              ex_not_found=.FALSE.
    149155           elseif((exlask(ilask).LT.min_ex ).and. ex_not_found ) then
    150               max_ex_iter=(min_ex-exlask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
    151                                / (exlask(ilask-1)-exlask(ilask))
     156              max_ex_iter=(min_ex-exlask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
     157                               / (exlask(ilask)-exlask(ilask-1))+yearlask(ilask-1)-Year
    152158              ex_not_found=.FALSE.
    153159           endif
     160           if(max_lsp_modulo .and. lsplask(ilask).lt.180) lsplask(ilask)=lsplask(ilask)+360.
     161           if(min_lsp_modulo .and. lsplask(ilask).gt.180) lsplask(ilask)=lsplask(ilask)-360.
    154162           if((lsplask(ilask).GT.max_lsp).and. lsp_not_found ) then
    155               max_lsp_iter=(max_lsp-lsplask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
    156                                / (lsplask(ilask-1)-lsplask(ilask))
     163              max_lsp_iter=(max_lsp-lsplask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
     164                               / (lsplask(ilask)-lsplask(ilask-1))+yearlask(ilask-1)-Year
    157165              lsp_not_found=.FALSE.
    158166           elseif((lsplask(ilask).LT.min_lsp ).and. lsp_not_found ) then
    159               max_lsp_iter=(min_lsp-lsplask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
    160                                / (lsplask(ilask-1)-lsplask(ilask))
     167              max_lsp_iter=(min_lsp-lsplask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) &
     168                               / (lsplask(ilask)-lsplask(ilask-1))+yearlask(ilask-1)-Year
    161169              lsp_not_found=.FALSE.
    162170           endif
     
    171179      print *, "Maximum number of iteration for the ex. parameter=", max_ex_iter
    172180
    173       print *, "Maximum lsp accepted=", max_obl
    174       print *, "Minimum lsp accepted=", min_obl
     181      print *, "Maximum lsp accepted=", max_lsp
     182      print *, "Minimum lsp accepted=", min_lsp
    175183      print *, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter
    176184
  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r2961 r2963  
    3232!module needed for INITIALISATION
    3333#ifndef CPP_STD
    34       use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
    35       use surfdat_h, only: tsurf, co2ice, emis,&
     34      use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa,inertiesoil
     35      use surfdat_h, only: tsurf, emis,&
    3636                           qsurf,watercap, ini_surfdat_h, &
    3737                           albedodat, zmea, zstd, zsig, zgam, zthe, &
     
    4040      use dimradmars_mod, only: totcloudfrac, albedo
    4141      use dust_param_mod, only: tauscaling
    42       use tracer_mod, only: noms,igcm_h2o_ice ! tracer names
     42      use tracer_mod, only: noms,igcm_h2o_ice,igcm_co2 ! tracer names
    4343#else
    4444      use comsoil_h, only: nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
    4545      use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe, &
    4646                           emissiv
    47       use tracer_h, only: noms,igcm_h2o_ice,igcm_co2_ice ! tracer names
     47      use tracer_h, only: noms,igcm_h2o_ice,igcm_co2 ! tracer names
    4848      use phys_state_var_mod
    4949#endif
     
    7373      USE mod_const_mpi, ONLY: COMM_LMDZ
    7474      USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, &
    75                            subslope_dist,co2ice_slope, &
    76                            tsurf_slope,tsoil_slope,fluxgrd_slope,&
    77                            fluxrad_sky_slope,sky_slope,callsubslope,&
    78                            co2iceflow, beta_slope, capcal_slope,&
    79                            albedo_slope,emiss_slope,qsurf_slope,&
    80                            iflat,major_slope,ini_comslope_h,fluxgeo_slope
     75                           subslope_dist,iflat,                &
     76                           major_slope,ini_comslope_h
    8177      use time_phylmdz_mod, only: daysec,dtphys
    8278      USE comconst_mod, ONLY: rad,g,r,cpp,pi
     
    9490                              TI_PEM,inertiedat_PEM,         &                        ! soil thermal inertia         
    9591                              tsoil_PEM, mlayer_PEM,layer_PEM,                  &     ! Soil temp, number of subsurface layers, soil mid layer depths
    96                               fluxgeo, &                             ! geothermal flux for the PEM and GCM
     92                              fluxgeo, &                                              ! geothermal flux for the PEM and GCM
    9793                              water_reservoir                                         ! Water ressources
    98       use adsorption_mod, only : regolith_adsorption,adsorption_pem, &   ! bool to check if adsorption, main subroutine
    99                                  ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! allocate arrays                               
    100                                  co2_adsorbded_phys, h2o_adsorbded_phys  ! mass of co2 and h2O adsorbded
     94      use adsorption_mod, only : regolith_adsorption,adsorption_pem, &                ! bool to check if adsorption, main subroutine
     95                                 ini_adsorption_h_PEM, end_adsorption_h_PEM, &        ! allocate arrays                               
     96                                 co2_adsorbded_phys, h2o_adsorbded_phys               ! mass of co2 and h2O adsorbded
    10197
    10298!!! For orbit parameters
     
    198194
    199195!!!!!!!!!!!!!!!!!!!!!!!! SLOPE
    200       REAL ,allocatable :: watercap_slope(:,:)                           ! Physics x Nslope: watercap per slope
    201       REAL ::              watercap_slope_saved                          ! Value saved at the previous time step
     196      REAL ::              watercap_saved                          ! Value saved at the previous time step
    202197      REAL , dimension(:,:), allocatable :: min_co2_ice_1                ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
    203198      REAL , dimension(:,:), allocatable :: min_co2_ice_2                ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
    204199      REAL , dimension(:,:), allocatable :: min_h2o_ice_1                ! ngrid field : minimum of water ice at each point for the first year [kg/m^2]
    205200      REAL , dimension(:,:), allocatable :: min_h2o_ice_2                ! ngrid field : minimum of water ice at each point for the second year [kg/m^2]
    206       REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_slope          ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
    207       REAL, dimension(:,:),allocatable  :: initial_co2_ice_sublim_slope    ! physical point field : Logical array indicating sublimating point of co2 ice
    208       REAL, dimension(:,:),allocatable  :: initial_h2o_ice_slope           ! physical point field : Logical array indicating if there is water ice at initial state
    209       REAL, dimension(:,:),allocatable  :: initial_co2_ice_slope           ! physical point field : Logical array indicating if there is co2 ice at initial state
     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
    210205      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice              ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
    211206      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
     
    225220     REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:)    !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    226221     REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
    227      REAL, ALLOCATABLE :: inertiesoil(:,:)    !Physic x Depth  Thermal inertia of the mesh for restart [SI]
    228      REAL, ALLOCATABLE :: TI_GCM(:,:,:) ! Same but for the start
     222
    229223     REAL,ALLOCATABLE  :: TI_locslope(:,:)    ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
    230224     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
     
    310304!----------------------------Initialisation : READ some constant of startfi_evol.nc ---------------------
    311305
    312 ! In the gcm, these values are given to the physic by the dynamic.
    313 ! Here we simply read them in the startfi_evol.nc file
    314       status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
    315 
    316       allocate(longitude(ngrid))
    317       allocate(latitude(ngrid))
    318       allocate(cell_area(ngrid))
    319 
    320       status = nf90_inq_varid(ncid, "longitude", lonvarid)
    321       status = nf90_get_var(ncid, lonvarid, longitude)
    322 
    323       status = nf90_inq_varid(ncid, "latitude", latvarid)
    324       status = nf90_get_var(ncid, latvarid, latitude)
    325 
    326       status = nf90_inq_varid(ncid, "area", areavarid)
    327       status = nf90_get_var(ncid, areavarid, cell_area)
    328 
    329       call ini_comsoil_h(ngrid)
    330 
    331       status = nf90_inq_varid(ncid, "soildepth", sdvarid)
    332       status = nf90_get_var(ncid, sdvarid, mlayer)
    333 
    334       status =nf90_close(ncid)
    335 
    336306!----------------------------READ start.nc ---------------------
    337307
     
    360330          iflag_phys)
    361331
     332! In the gcm, these values are given to the physic by the dynamic.
     333! Here we simply read them in the startfi_evol.nc file
     334      status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
     335
     336      allocate(longitude(ngrid))
     337      allocate(latitude(ngrid))
     338      allocate(cell_area(ngrid))
     339
     340      status = nf90_inq_varid(ncid, "longitude", lonvarid)
     341      status = nf90_get_var(ncid, lonvarid, longitude)
     342
     343      status = nf90_inq_varid(ncid, "latitude", latvarid)
     344      status = nf90_get_var(ncid, latvarid, latitude)
     345
     346      status = nf90_inq_varid(ncid, "area", areavarid)
     347      status = nf90_get_var(ncid, areavarid, cell_area)
     348
     349      status = nf90_inq_varid(ncid, "soildepth", sdvarid)
     350      status = nf90_get_var(ncid, sdvarid, mlayer)
     351
     352      status =nf90_close(ncid)
     353
    362354!----------------------------READ startfi.nc ---------------------
    363355
    364356! First we read the initial state (starfi.nc)
    365 
    366       allocate(watercap_slope(ngrid,nslope))
    367       allocate(TI_GCM(ngrid,nsoilmx,nslope))
    368       allocate(inertiesoil(ngrid,nsoilmx))
    369357
    370358#ifndef CPP_STD
     
    373361              day_ini,time_phys,         &
    374362              tsurf,tsoil,albedo,emis,   &
    375               q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar,     &
    376               watercap,inertiesoil,nslope,tsurf_slope,           &
    377               tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
    378               subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM,     &
    379               qsurf_slope,watercap_slope)
    380 
    381 
     363              q2,qsurf,tauscaling,totcloudfrac,wstar,     &
     364              watercap,def_slope,def_slope_mean,subslope_dist)
    382365
    383366! Remove unphysical values of surface tracer
     
    385368       DO nnq=1,nqtot
    386369         DO islope=1,nslope
    387            if(qsurf_slope(i,nnq,islope).LT.0) then
    388              qsurf_slope(i,nnq,islope)=0.
     370           if(qsurf(i,nnq,islope).LT.0) then
     371             qsurf(i,nnq,islope)=0.
    389372           endif
    390373         enddo
     
    409392
    410393         allocate(co2ice(ngrid))
    411          co2ice(:)=qsurf(:,igcm_co2_ice)
     394         co2ice(:)=qsurf(:,igcm_co2)
    412395         co2ice_slope(:,1)=co2ice(:)
    413396         tsurf_slope(:,1)=tsurf(:)
     
    433416     DO nnq=1,nqtot
    434417       if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
     418       if(noms(nnq).eq."co2") igcm_co2 = nnq
    435419     ENDDO
    436420
     
    492476     allocate(q_co2_PEM_phys(ngrid,timelen))
    493477     allocate(q_h2o_PEM_phys(ngrid,timelen))
    494      allocate(co2_ice_GCM_slope(ngrid,nslope,timelen))
     478     allocate(co2_ice_GCM(ngrid,nslope,timelen))
    495479     allocate(watersurf_density_ave(ngrid,nslope))
    496480     allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
     
    508492
    509493     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,&   
    510                        tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM_slope, &     
     494                       tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, &     
    511495                       watersurf_density_ave,watersoil_density_timeseries)
    512496
     
    517501
    518502     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, &
    519                   tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM_slope, &     
     503                  tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, &     
    520504                  watersurf_density_ave,watersoil_density_timeseries)
    521505
     
    543527
    544528   if(soil_pem) then
    545       call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM,TI_PEM)
     529      call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
    546530      DO l=1,nsoilmx
    547531        tsoil_ave_PEM_yr1(:,l,:)=tsoil_ave_yr1(:,l,:)
     
    605589!---------------------------- Save initial PCM situation ---------------------
    606590
    607      allocate(initial_co2_ice_sublim_slope(ngrid,nslope))
    608      allocate(initial_co2_ice_slope(ngrid,nslope))
    609      allocate(initial_h2o_ice_slope(ngrid,nslope))
     591     allocate(initial_co2_ice_sublim(ngrid,nslope))
     592     allocate(initial_co2_ice(ngrid,nslope))
     593     allocate(initial_h2o_ice(ngrid,nslope))
    610594
    611595! We save the places where water ice is sublimating
     
    618602    do islope=1,nslope
    619603      if (tendencies_co2_ice(i,islope).LT.0) then
    620          initial_co2_ice_sublim_slope(i,islope)=1.
     604         initial_co2_ice_sublim(i,islope)=1.
    621605         ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
    622606      else
    623          initial_co2_ice_sublim_slope(i,islope)=0.         
     607         initial_co2_ice_sublim(i,islope)=0.         
    624608      endif
    625       if (co2ice_slope(i,islope).GT.0) then
    626          initial_co2_ice_slope(i,islope)=1.
     609      if (qsurf(i,igcm_co2,islope).GT.0) then
     610         initial_co2_ice(i,islope)=1.
    627611      else
    628          initial_co2_ice_slope(i,islope)=0.         
     612         initial_co2_ice(i,islope)=0.         
    629613      endif
    630614      if (tendencies_h2o_ice(i,islope).LT.0) then
    631          initial_h2o_ice_slope(i,islope)=1.
     615         initial_h2o_ice(i,islope)=1.
    632616         ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
    633617      else
    634          initial_h2o_ice_slope(i,islope)=0.         
     618         initial_h2o_ice(i,islope)=0.         
    635619      endif
    636620    enddo
     
    674658      call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_PEM_yr1,tsoil_PEM,porefillingice_depth,porefillingice_thickness, &
    675659      tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,&
    676       tendencies_h2o_ice,tendencies_co2_ice,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,&
     660      tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,&
    677661      watersurf_density_ave,watersoil_density_PEM_ave, &
    678662      co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
     
    680664      do ig = 1,ngrid
    681665        do islope = 1,nslope
    682           qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
     666          qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    683667        enddo
    684668      enddo
     
    829813      print *, "Evolution of h2o ice"
    830814     
    831      call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
     815     call evol_h2o_ice_s_slope(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
    832816     
    833817     DO islope=1, nslope
    834818       write(str2(1:2),'(i2.2)') islope
    835        call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf_slope(:,igcm_h2o_ice,islope))
     819       call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
    836820       call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
    837821     ENDDO
    838822
    839823      print *, "Evolution of co2 ice"
    840      call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
     824     call evol_co2_ice_s_slope(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
    841825
    842826!------------------------
     
    853837      if (co2glaciersflow) then
    854838       call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,&
    855                          global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
     839                         global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
    856840      endif
    857841
     
    860844     DO islope=1, nslope
    861845       write(str2(1:2),'(i2.2)') islope
    862        call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,co2ice_slope(:,islope))
     846       call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
    863847       call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
    864848     ENDDO
     
    881865      DO ig = 1,ngrid
    882866        DO islope = 1,nslope
    883           if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-10) THEN !co2ice disappeared, look for closest point without co2ice
     867          if(initial_co2_ice(ig,islope).gt.0.5 .and. qsurf(ig,igcm_co2,islope).LT. 1E-10) THEN !co2ice disappeared, look for closest point without co2ice
    884868            if(latitude_deg(ig).gt.0) then
    885869              do ig_loop=ig,ngrid
    886870                DO islope_loop=islope,iflat,-1
    887                   if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
     871                  if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then
    888872                    tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
    889873                    bool_sublim=1
     
    898882              do ig_loop=ig,1,-1
    899883                DO islope_loop=islope,iflat
    900                   if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
     884                  if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then
    901885                    tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
    902886                    bool_sublim=1
     
    909893              enddo
    910894            endif
    911             initial_co2_ice_slope(ig,islope)=0
    912             if ((co2ice_slope(ig,islope).lt.1e-10).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then   
    913               albedo_slope(ig,1,islope) = albedo_h2o_frost
    914               albedo_slope(ig,2,islope) = albedo_h2o_frost
     895            initial_co2_ice(ig,islope)=0
     896            if ((qsurf(ig,igcm_co2,islope).lt.1e-10).and. (qsurf(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then   
     897              albedo(ig,1,islope) = albedo_h2o_frost
     898              albedo(ig,2,islope) = albedo_h2o_frost
    915899            else
    916               albedo_slope(ig,1,islope) = albedodat(ig)
    917               albedo_slope(ig,2,islope) = albedodat(ig)
    918               emiss_slope(ig,islope) = emissiv         
     900              albedo(ig,1,islope) = albedodat(ig)
     901              albedo(ig,2,islope) = albedodat(ig)
     902              emis(ig,islope) = emissiv         
    919903            endif
    920           elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
     904          elseif( (qsurf(ig,igcm_co2,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
    921905            ave=0.
    922906            do t=1,timelen
    923               if(co2_ice_GCM_slope(ig,islope,t).gt.1e-3) then
     907              if(co2_ice_GCM(ig,islope,t).gt.1e-3) then
    924908                ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
    925909              else
     
    938922      do ig = 1,ngrid
    939923        do islope = 1,nslope
    940           tsurf_slope(ig,islope) =  tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope))
     924          tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope))
    941925        enddo
    942926      enddo
     
    944928     DO islope=1, nslope
    945929       write(str2(1:2),'(i2.2)') islope
    946        call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_slope(:,islope))
     930       call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
    947931     ENDDO
    948932
     
    990974
    991975! II_d.4 Update the soil thermal properties
    992       call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
     976      call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, &
    993977        porefillingice_depth,porefillingice_thickness,TI_PEM)
    994978
    995979! II_d.5 Update the mass of the regolith adsorbded
    996980     if(adsorption_pem) then
    997        call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf_slope(:,igcm_h2o_ice,:),co2ice_slope, &
     981       call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:), &
    998982                             tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
    999983                             h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
     
    1013997
    1014998      print *, "Adaptation of the new co2 tendencies to the current pressure"
    1015       call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,co2ice_slope,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,&     
     999      call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,&     
    10161000                               global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
    10171001
     
    10271011
    10281012!------------------------
    1029       call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
    1030 
    1031       call criterion_co2_stop(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,STOPPING_pressure,ngrid, &
    1032                                    initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
     1013      call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
     1014
     1015      call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
     1016                                   initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
    10331017
    10341018      year_iter=year_iter+dt_pem     
     
    10841068
    10851069! III_a.1 Ice update (for startfi)
     1070#ifdef CPP_STD
    10861071! Co2 ice
    10871072      DO ig = 1,ngrid
     
    10921077                      cos(pi*def_slope_mean(islope)/180.)
    10931078        ENDDO
    1094 #ifdef CPP_STD
    1095         qsurf(ig,igcm_co2_ice)=co2ice(ig)
     1079        qsurf(ig,igcm_co2)=co2ice(ig)
     1080      ENDDO ! of DO ig=1,ngrid
    10961081#endif
    1097       ENDDO ! of DO ig=1,ngrid
    10981082
    10991083! H2O ice
     
    11021086         watercap_sum=0.
    11031087         DO islope=1,nslope
    1104            watercap_slope_saved = watercap_slope(ig,islope)
    1105            if(qsurf_slope(ig,igcm_h2o_ice,islope).GT. (watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then
    1106              qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
     1088           watercap_saved = watercap(ig,islope)
     1089           if(qsurf(ig,igcm_h2o_ice,islope).GT. (watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then
     1090             qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
    11071091           else
    1108              watercap_slope(ig,islope)=watercap_slope(ig,islope)+qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
    1109              qsurf_slope(ig,igcm_h2o_ice,islope)=0.
     1092             watercap(ig,islope)=watercap(ig,islope)+qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
     1093             qsurf(ig,igcm_h2o_ice,islope)=0.
    11101094           endif
    1111            watercap_sum=watercap_sum+(watercap_slope(ig,islope)-watercap_slope_saved)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    1112            watercap_slope(ig,islope)=0.
     1095           watercap_sum=watercap_sum+(watercap(ig,islope)-watercap_saved)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
     1096           watercap(ig,islope)=0.
    11131097         enddo
    11141098           water_reservoir(ig)=water_reservoir(ig)+watercap_sum
     
    11161100     enddo
    11171101
    1118       DO ig = 1,ngrid
    1119         qsurf(ig,igcm_h2o_ice) = 0.
    1120         DO islope = 1,nslope
    1121           qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
    1122                                    * subslope_dist(ig,islope) /          &
    1123                                   cos(pi*def_slope_mean(islope)/180.)
    1124         ENDDO
    1125       ENDDO ! of DO ig=1,ngrid
    1126 
    11271102     DO ig=1,ngrid
    11281103       DO islope=1,nslope
    1129          if(qsurf_slope(ig,igcm_h2o_ice,islope).GT.500) then
     1104         if(qsurf(ig,igcm_h2o_ice,islope).GT.500) then
    11301105            watercaptag(ig)=.true.
    1131             qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-250
     1106            qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)-250
    11321107            water_reservoir(ig)=water_reservoir(ig)+250
    11331108         endif
     
    11361111
    11371112     DO ig=1,ngrid
    1138          if(water_reservoir(ig).LT. 10) then
    1139             watercaptag(ig)=.false.
    1140             qsurf(ig,igcm_h2o_ice)=qsurf(ig,igcm_h2o_ice)+water_reservoir(ig)
    1141        DO islope=1,nslope
    1142             qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+water_reservoir(ig)
    1143        ENDDO
    1144          endif
    1145      enddo
    1146 
    1147       DO ig = 1,ngrid
    1148         watercap(ig) = 0.
    1149         DO islope = 1,nslope
    1150           watercap(ig) = watercap(ig) + watercap_slope(ig,islope) &
    1151                                    * subslope_dist(ig,islope) /          &
    1152                                   cos(pi*def_slope_mean(islope)/180.)
    1153         ENDDO
    1154       ENDDO ! of DO ig=1,ngrid
     1113       if(water_reservoir(ig).LT. 10) then
     1114         watercaptag(ig)=.false.
     1115         DO islope=1,nslope
     1116           qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+water_reservoir(ig)
     1117         ENDDO
     1118       endif
     1119     enddo
    11551120
    11561121! III_a.2  Tsoil update (for startfi)
    11571122
    11581123   if(soil_pem) then
    1159      fluxgeo_slope(:,:) = fluxgeo
    1160      call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM)
    1161      tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)     
     1124     call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
     1125     tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)     
    11621126   endif !soil_pem
    1163 
    1164 #ifndef CPP_STD
    1165       DO ig = 1,ngrid
    1166         DO iloop = 1,nsoilmx
    1167           tsoil(ig,iloop) = 0.
    1168           inertiesoil(ig,iloop) = 0.
    1169           DO islope = 1,nslope
    1170             tsoil(ig,iloop) =  tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
    1171                             * subslope_dist(ig,islope)   
    1172             inertiesoil(ig,iloop) =  inertiesoil(ig,iloop) + TI_GCM(ig,iloop,islope) &
    1173                             * subslope_dist(ig,islope)         
    1174           ENDDO
    1175         ENDDO
    1176       ENDDO ! of DO ig=1,ngrid
    1177 
    1178 ! III_a.3 Surface optical properties (for startfi)
    1179 
    1180     DO ig = 1,ngrid
    1181        DO l = 1,2
    1182         albedo(ig,l) =0.
    1183         DO islope = 1,nslope
    1184           albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) &
    1185                          *subslope_dist(ig,islope)
    1186         ENDDO
    1187       ENDDO
    1188     ENDDO
    1189 
    1190     DO ig = 1,ngrid
    1191       emis(ig) =0.
    1192       DO islope = 1,nslope
    1193         emis(ig)= emis(ig)+emiss_slope(ig,islope) &
    1194                          *subslope_dist(ig,islope)
    1195       ENDDO
    1196     ENDDO
    1197 
    1198       DO ig = 1,ngrid
    1199         tsurf(ig) = 0.
    1200         DO islope = 1,nslope
    1201           tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)*tsurf_slope(ig,islope))**4 &
    1202                      * subslope_dist(ig,islope)     
    1203         ENDDO
    1204         tsurf(ig) = tsurf(ig)**(0.25)/emis(ig)
    1205       ENDDO ! of DO ig=1,ngrid
    1206 
    1207 #endif
    12081127
    12091128! III_a.4 Pressure (for start)
     
    12501169       DO ig=1,ngrid
    12511170         DO l=1,llm-1
    1252            if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") ) then
     1171           if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") .and. (noms(nnq).NE."stormdust_number") .and. (noms(nnq).NE."topdust_number")) then
    12531172             extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
    12541173             q(ig,l,nnq)=1.
     
    12951214                        nsoilmx,ngrid,nlayer,nq,              &
    12961215                        ptimestep,pday,0.,cell_area,          &
    1297                         albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
    1298                         hmons,summit,base,nslope,def_slope,   &
     1216                        albedodat,inertiedat,def_slope,   &
    12991217                        subslope_dist)
    13001218
    13011219     call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
    13021220                     ptimestep,ztime_fin,                        &
    1303                      tsurf,tsoil,co2ice,albedo,emis,             &
    1304                      q2,qsurf,tauscaling,totcloudfrac,wstar,     &
    1305                      watercap,inertiesoil,nslope,co2ice_slope,   &
    1306                      tsurf_slope,tsoil_slope, albedo_slope,      &
    1307                      emiss_slope,qsurf_slope,watercap_slope, TI_GCM)
     1221                     tsurf,tsoil,inertiesoil,albedo,             &
     1222                     emis,q2,qsurf,tauscaling,totcloudfrac,wstar,&
     1223                     watercap)
    13081224#else
    13091225     call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
     
    13461262     deallocate(q_co2_PEM_phys)
    13471263     deallocate(q_h2o_PEM_phys)
    1348      deallocate(co2_ice_GCM_slope)
     1264     deallocate(co2_ice_GCM)
    13491265     deallocate(watersurf_density_ave)
    13501266     deallocate(watersoil_density_timeseries)
  • trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90

    r2944 r2963  
    9191     print *, "Downloading data for vmr co2..."
    9292
    93   CALL get_var3("co2_cropped"   ,q_co2_dyn)
     93  CALL get_var3("co2_layer1"   ,q_co2_dyn)
    9494
    9595     print *, "Downloading data for vmr co2 done"
    9696     print *, "Downloading data for vmr h20..."
    9797
    98   CALL get_var3("h2o_cropped"   ,q_h2o_dyn)
     98  CALL get_var3("h2o_layer1"   ,q_h2o_dyn)
    9999
    100100     print *, "Downloading data for vmr h2o done"
     
    143143     if(soil_pem) then
    144144
    145      print *, "Downloading data for tsoil_slope ..."
    146 
    147 DO islope=1,nslope
    148   write(num,fmt='(i2.2)') islope
    149   call get_var4("tsoil_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
    150 ENDDO
    151 
    152      print *, "Downloading data for tsoil_slope done"
     145     print *, "Downloading data for soiltemp_slope ..."
     146
     147DO islope=1,nslope
     148  write(num,fmt='(i2.2)') islope
     149  call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
     150ENDDO
     151
     152     print *, "Downloading data for soiltemp_slope done"
    153153
    154154     print *, "Downloading data for watersoil_density ..."
     
    182182
    183183    if(soil_pem) then
    184       call get_var4("tsoil",tsoil_gcm_dyn(:,:,:,1,:))
     184      call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
    185185    endif !soil_pem
    186186  endif !nslope=1
  • trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90

    r2894 r2963  
    1616#ifndef CPP_STD     
    1717      USE comconst_mod, ONLY: pi
    18       USE planete_h, ONLY: e_elips, obliquit, timeperi
     18      USE planete_h, ONLY: e_elips, obliquit, timeperi, periheli,aphelie,p_elips,peri_day,year_day
    1919#else
    2020      use planete_mod, only: e_elips, obliquit, timeperi
     
    4545      real :: timeperi_ls               ! time of peri in ls
    4646      integer ilask                     ! Loop variable
     47      real :: halfaxe                   ! Million of km
     48      real :: unitastr
    4749
    4850
     
    7375             endif
    7476             if(var_lsp) then
     77               if(lsplask(ilask)-lsplask(ilask+1) .gt.200) then
     78                 if(lsplask(ilask).lt.lsplask(ilask+1)) then
     79                   lsplask(ilask)=lsplask(ilask)+360
     80                 else
     81                   lsplask(ilask+1)=lsplask(ilask+1)+360
     82                 endif
     83               endif
    7584               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.)
    7686             endif
    7787              exit
    7888           endif
    7989        enddo
     90        halfaxe=227.94
    8091        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)
    8197
    8298       print *, "recomp_orb_param, Final year of the PEM run:", Year
  • trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90

    r2855 r2963  
    5151  integer :: numyear
    5252
    53 
    54 
    55 
    56 
    5753DO numyear=1, 2
    5854write(*,*) 'numyear',numyear
     
    6157!integer :: ncid, status
    6258!...
    63 status = nf90_open(path = "Xdiurnalave"//str2//".nc", mode = nf90_nowrite, ncid = ncid1)
    64 if(status /= nf90_noerr) call handle_err(status)
    65 
    66 status = nf90_create(path = "Xdiurnalave"//str2//".nc_new", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2)
     59status = nf90_open(path = "data2reshape"//str2//".nc", mode = nf90_nowrite, ncid = ncid1)
     60if(status /= nf90_noerr) call handle_err(status)
     61
     62status = nf90_create(path = "datareshaped"//str2//".nc", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2)
    6763if(status /= nf90_noerr) call handle_err(status)
    6864
  • trunk/LMDZ.MARS/deftank/context_lmdz_physics.xml

    r2934 r2963  
    129129
    130130        <!-- diurnal averages outputs; output_freq is every day -->
    131         <file id="diurnalave"
    132               name="Xdiurnalave"
     131        <file id="diurnalave_"
     132              name="Xdiurnalave_"
    133133              output_freq="1d"
    134134              time_units="days"
     
    157157            </field_group>
    158158        </file>
     159
     160        <file id="diurnalave"
     161              name="Xdiurnalave"
     162              output_freq="1d"
     163              type="one_file"
     164              time_units="days"
     165              enabled=".false.">
     166                   
     167            <!-- VARS 0D -->
     168            <field_group operation="average"
     169                         freq_op="1ts">
     170                <field field_ref="Ls" />
     171            </field_group>
     172
     173            <field_group operation="average"
     174                         freq_op="1ts">
     175                <field field_ref="area" operation="once" />
     176                <field field_ref="ps" />
     177                <field field_ref="tsurf" />
     178                <field field_ref="tsurf_slope01" />
     179                <field field_ref="tsurf_slope02" />
     180                <field field_ref="tsurf_slope03" />
     181                <field field_ref="tsurf_slope04" />
     182                <field field_ref="tsurf_slope05" />
     183                <field field_ref="tsurf_slope06" />
     184                <field field_ref="tsurf_slope07" />
     185                <field field_ref="Waterdensity_surface01" />
     186                <field field_ref="Waterdensity_surface02" />
     187                <field field_ref="Waterdensity_surface03" />
     188                <field field_ref="Waterdensity_surface04" />
     189                <field field_ref="Waterdensity_surface05" />
     190                <field field_ref="Waterdensity_surface06" />
     191                <field field_ref="Waterdensity_surface07" />
     192                <field field_ref="h2o_layer1" />
     193                <field field_ref="co2_layer1" />
     194            </field_group>
     195
     196            <field_group operation="minimum"
     197                         freq_op="1ts">
     198                <field field_ref="co2ice" />
     199                <field field_ref="co2ice_slope01" />
     200                <field field_ref="co2ice_slope02" />
     201                <field field_ref="co2ice_slope03" />
     202                <field field_ref="co2ice_slope04" />
     203                <field field_ref="co2ice_slope05" />
     204                <field field_ref="co2ice_slope06" />
     205                <field field_ref="co2ice_slope07" />
     206                <field field_ref="h2o_ice_s" />
     207                <field field_ref="h2o_ice_s_slope01" />
     208                <field field_ref="h2o_ice_s_slope02" />
     209                <field field_ref="h2o_ice_s_slope03" />
     210                <field field_ref="h2o_ice_s_slope04" />
     211                <field field_ref="h2o_ice_s_slope05" />
     212                <field field_ref="h2o_ice_s_slope06" />
     213                <field field_ref="h2o_ice_s_slope07" />
     214                <field field_ref="watercap_slope01" />
     215                <field field_ref="watercap_slope02" />
     216                <field field_ref="watercap_slope03" />
     217                <field field_ref="watercap_slope04" />
     218                <field field_ref="watercap_slope05" />
     219                <field field_ref="watercap_slope06" />
     220                <field field_ref="watercap_slope07" />
     221            </field_group>
     222
     223            <!-- VARS soil -->
     224            <field_group operation="average"
     225                         freq_op="1ts">
     226                <field field_ref="soiltemp" />
     227                <field field_ref="soiltemp_slope01" />
     228                <field field_ref="soiltemp_slope02" />
     229                <field field_ref="soiltemp_slope03" />
     230                <field field_ref="soiltemp_slope04" />
     231                <field field_ref="soiltemp_slope05" />
     232                <field field_ref="soiltemp_slope06" />
     233                <field field_ref="soiltemp_slope07" />
     234                <field field_ref="Waterdensity_soil_slope01" />
     235                <field field_ref="Waterdensity_soil_slope02" />
     236                <field field_ref="Waterdensity_soil_slope03" />
     237                <field field_ref="Waterdensity_soil_slope04" />
     238                <field field_ref="Waterdensity_soil_slope05" />
     239                <field field_ref="Waterdensity_soil_slope06" />
     240                <field field_ref="Waterdensity_soil_slope07" />
     241            </field_group>
     242        </file>
     243
     244
    159245    </file_definition>
    160246
  • trunk/LMDZ.MARS/deftank/field_def_physics_mars.xml

    r2934 r2963  
    6868                   long_name="Surface emissivity of slope 01" 
    6969                   unit="" />
     70            <field id="emis_slope02"
     71                   long_name="Surface emissivity of slope 02" 
     72                   unit="" />
     73            <field id="emis_slope03"
     74                   long_name="Surface emissivity of slope 03" 
     75                   unit="" />
     76            <field id="emis_slope04"
     77                   long_name="Surface emissivity of slope 04" 
     78                   unit="" />
     79            <field id="emis_slope05"
     80                   long_name="Surface emissivity of slope 05" 
     81                   unit="" />
     82            <field id="emis_slope06"
     83                   long_name="Surface emissivity of slope 06" 
     84                   unit="" />
     85            <field id="emis_slope07"
     86                   long_name="Surface emissivity of slope 07" 
     87                   unit="" />
    7088            <field id="albedo"
    7189                   long_name="Albedo of the surface" 
     
    7492                   long_name="Albedo of the surface for slope 01" 
    7593                   unit="" />
    76 
     94            <field id="albedo_slope02"
     95                   long_name="Albedo of the surface for slope 02" 
     96                   unit="" />
     97            <field id="albedo_slope03"
     98                   long_name="Albedo of the surface for slope 03" 
     99                   unit="" />
     100            <field id="albedo_slope04"
     101                   long_name="Albedo of the surface for slope 04" 
     102                   unit="" />
     103            <field id="albedo_slope05"
     104                   long_name="Albedo of the surface for slope 05" 
     105                   unit="" />
     106            <field id="albedo_slope06"
     107                   long_name="Albedo of the surface for slope 06" 
     108                   unit="" />
     109            <field id="albedo_slope07"
     110                   long_name="Albedo of the surface for slope 07" 
     111                   unit="" />
    77112            <field id="local_time"
    78113                   long_name="Local time"
     
    88123                   long_name="Surface Temperature of slope 01"
    89124                   unit="K" />
     125            <field id="tsurf_slope02"
     126                   long_name="Surface Temperature of slope 02"
     127                   unit="K" />
     128            <field id="tsurf_slope03"
     129                   long_name="Surface Temperature of slope 03"
     130                   unit="K" />
     131            <field id="tsurf_slope04"
     132                   long_name="Surface Temperature of slope 04"
     133                   unit="K" />
     134            <field id="tsurf_slope05"
     135                   long_name="Surface Temperature of slope 05"
     136                   unit="K" />
     137            <field id="tsurf_slope06"
     138                   long_name="Surface Temperature of slope 06"
     139                   unit="K" />
     140            <field id="tsurf_slope07"
     141                   long_name="Surface Temperature of slope 07"
     142                   unit="K" />
    90143            <field id="temp_layer1"
    91144                   long_name="Temperature in layer 1"
     
    102155                   long_name="Longwave radiation at the surface on slope 01"
    103156                   unit="W.m-2" />
     157            <field id="fluxsurf_lw_slope02"
     158                   long_name="Longwave radiation at the surface on slope 02"
     159                   unit="W.m-2" />
     160            <field id="fluxsurf_lw_slope03"
     161                   long_name="Longwave radiation at the surface on slope 03"
     162                   unit="W.m-2" />
     163            <field id="fluxsurf_lw_slope04"
     164                   long_name="Longwave radiation at the surface on slope 04"
     165                   unit="W.m-2" />
     166            <field id="fluxsurf_lw_slope05"
     167                   long_name="Longwave radiation at the surface on slope 05"
     168                   unit="W.m-2" />
     169            <field id="fluxsurf_lw_slope06"
     170                   long_name="Longwave radiation at the surface on slope 06"
     171                   unit="W.m-2" />
     172            <field id="fluxsurf_lw_slope07"
     173                   long_name="Longwave radiation at the surface on slope 07"
     174                   unit="W.m-2" />
    104175            <field id="fluxtop_lw"
    105176                   long_name="Longwave radiation at the top of the atmosphere"
     
    111182            <field id="fluxsurf_dn_sw_slope01"
    112183                   long_name="Incoming shortwave radiation at the surface on slope 01" 
     184                   unit="W.m-2" />
     185            <field id="fluxsurf_dn_sw_slope02"
     186                   long_name="Incoming shortwave radiation at the surface on slope 02" 
     187                   unit="W.m-2" />
     188            <field id="fluxsurf_dn_sw_slope03"
     189                   long_name="Incoming shortwave radiation at the surface on slope 03" 
     190                   unit="W.m-2" />
     191            <field id="fluxsurf_dn_sw_slope04"
     192                   long_name="Incoming shortwave radiation at the surface on slope 04" 
     193                   unit="W.m-2" />
     194            <field id="fluxsurf_dn_sw_slope05"
     195                   long_name="Incoming shortwave radiation at the surface on slope 05" 
     196                   unit="W.m-2" />
     197            <field id="fluxsurf_dn_sw_slope06"
     198                   long_name="Incoming shortwave radiation at the surface on slope 06" 
     199                   unit="W.m-2" />
     200            <field id="fluxsurf_dn_sw_slope07"
     201                   long_name="Incoming shortwave radiation at the surface on slope 07" 
    113202                   unit="W.m-2" />
    114203            <field id="fluxtop_up_sw"
     
    145234                   long_name="tendency due to vertical diffusion of background dust on surface"
    146235                   unit="kg.m-2.s-1" />
     236            <field id="dqsdifdustq_slope02"
     237                   long_name="tendency due to vertical diffusion of background dust on surface"
     238                   unit="kg.m-2.s-1" />
     239            <field id="dqsdifdustq_slope03"
     240                   long_name="tendency due to vertical diffusion of background dust on surface"
     241                   unit="kg.m-2.s-1" />
     242            <field id="dqsdifdustq_slope04"
     243                   long_name="tendency due to vertical diffusion of background dust on surface"
     244                   unit="kg.m-2.s-1" />
     245            <field id="dqsdifdustq_slope05"
     246                   long_name="tendency due to vertical diffusion of background dust on surface"
     247                   unit="kg.m-2.s-1" />
     248            <field id="dqsdifdustq_slope06"
     249                   long_name="tendency due to vertical diffusion of background dust on surface"
     250                   unit="kg.m-2.s-1" />
     251            <field id="dqsdifdustq_slope07"
     252                   long_name="tendency due to vertical diffusion of background dust on surface"
     253                   unit="kg.m-2.s-1" />
    147254            <field id="dqsdifrdsq"
    148255                   long_name="tendency due to vertical diffusion of stormdust on surface"
    149256                   unit="kg.m-2.s-1" />
    150257            <field id="dqsdifrdsq_slope01"
     258                   long_name="tendency due to vertical diffusion of stormdust on surface"
     259                   unit="kg.m-2.s-1" />
     260            <field id="dqsdifrdsq_slope02"
     261                   long_name="tendency due to vertical diffusion of stormdust on surface"
     262                   unit="kg.m-2.s-1" />
     263            <field id="dqsdifrdsq_slope03"
     264                   long_name="tendency due to vertical diffusion of stormdust on surface"
     265                   unit="kg.m-2.s-1" />
     266            <field id="dqsdifrdsq_slope04"
     267                   long_name="tendency due to vertical diffusion of stormdust on surface"
     268                   unit="kg.m-2.s-1" />
     269            <field id="dqsdifrdsq_slope05"
     270                   long_name="tendency due to vertical diffusion of stormdust on surface"
     271                   unit="kg.m-2.s-1" />
     272            <field id="dqsdifrdsq_slope06"
     273                   long_name="tendency due to vertical diffusion of stormdust on surface"
     274                   unit="kg.m-2.s-1" />
     275            <field id="dqsdifrdsq_slope07"
    151276                   long_name="tendency due to vertical diffusion of stormdust on surface"
    152277                   unit="kg.m-2.s-1" />
     
    185310            <field id="watercap_slope01"
    186311                   long_name="Perennial water ice thickness of slope 01"
    187                    unit="kg.m-2" />           
     312                   unit="kg.m-2" /> 
     313            <field id="watercap_slope02"
     314                   long_name="Perennial water ice thickness of slope 02"
     315                   unit="kg/m2" />
     316            <field id="watercap_slope03"
     317                   long_name="Perennial water ice thickness of slope 03"
     318                   unit="kg/m2" />
     319            <field id="watercap_slope04"
     320                   long_name="Perennial water ice thickness of slope 04"
     321                   unit="kg/m2" />
     322            <field id="watercap_slope05"
     323                   long_name="Perennial water ice thickness of slope 05"
     324                   unit="kg/m2" />
     325            <field id="watercap_slope06"
     326                   long_name="Perennial water ice thickness of slope 06"
     327                   unit="kg/m2" />
     328            <field id="watercap_slope07"
     329                   long_name="Perennial water ice thickness of slope 07"
     330                   unit="kg/m2" />
     331   
    188332            <field id="surf_h2o_lh"
    189333                   long_name="Ground ice latent heat flux"
     
    213357                   long_name="Mass of water ice on the surface of slope 01"
    214358                   unit="kg.m-2" />
     359            <field id="h2o_ice_s_slope02"
     360                   long_name="Mass of water ice on the surface of slope 02"
     361                   unit="kg/m2" />
     362            <field id="h2o_ice_s_slope03"
     363                   long_name="Mass of water ice on the surface of slope 03"
     364                   unit="kg/m2" />
     365            <field id="h2o_ice_s_slope04"
     366                   long_name="Mass of water ice on the surface of slope 04"
     367                   unit="kg/m2" />
     368            <field id="h2o_ice_s_slope05"
     369                   long_name="Mass of water ice on the surface of slope 05"
     370                   unit="kg/m2" />
     371            <field id="h2o_ice_s_slope06"
     372                   long_name="Mass of water ice on the surface of slope 06"
     373                   unit="kg/m2" />
     374            <field id="h2o_ice_s_slope07"
     375                   long_name="Mass of water ice on the surface of slope 07"
     376                   unit="kg/m2" />
     377            <field id="Waterdensity_surface01"
     378                   long_name="Waterdensity_surface of slope 01"
     379                   unit="XX" />
     380            <field id="Waterdensity_surface02"
     381                   long_name="Waterdensity_surface of slope 02"
     382                   unit="XX" />
     383            <field id="Waterdensity_surface03"
     384                   long_name="Waterdensity_surface of slope 03"
     385                   unit="XX" />
     386            <field id="Waterdensity_surface04"
     387                   long_name="Waterdensity_surface of slope 04"
     388                   unit="XX" />
     389            <field id="Waterdensity_surface05"
     390                   long_name="Waterdensity_surface of slope 05"
     391                   unit="XX" />
     392            <field id="Waterdensity_surface06"
     393                   long_name="Waterdensity_surface of slope 06"
     394                   unit="XX" />
     395            <field id="Waterdensity_surface07"
     396                   long_name="Waterdensity_surface of slope 07"
     397                   unit="XX" />
     398            <field id="h2o_layer1"
     399                   long_name="h2o in the first layer" 
     400                   unit="kg/kg" />
     401            <field id="co2_layer1"
     402                   long_name="co2 in the first layer" 
     403                   unit="kg/kg" />
     404
     405
    215406           
    216407            <!-- CO2 cycle -->
     
    221412                   long_name="CO2 ice thickness on slope 01"
    222413                   unit="kg.m-2" />
     414            <field id="co2ice_slope02"
     415                   long_name="CO2 ice thickness on slope 02"
     416                   unit="kg/m2" />
     417            <field id="co2ice_slope03"
     418                   long_name="CO2 ice thickness on slope 03"
     419                   unit="kg/m2" />
     420            <field id="co2ice_slope04"
     421                   long_name="CO2 ice thickness on slope 04"
     422                   unit="kg/m2" />
     423            <field id="co2ice_slope05"
     424                   long_name="CO2 ice thickness on slope 05"
     425                   unit="kg/m2" />
     426            <field id="co2ice_slope06"
     427                   long_name="CO2 ice thickness on slope 06"
     428                   unit="kg/m2" />
     429            <field id="co2ice_slope07"
     430                   long_name="CO2 ice thickness on slope 07"
     431                   unit="kg/m2" />
     432
    223433           
    224434            <field id="co2condens_zfallice"
     
    463673                   long_name="Soil temperature for slope 01"
    464674                   unit="K" />
     675            <field id="soiltemp_slope02"
     676                   long_name="Soil temperature for slope 02"
     677                   unit="K" />
     678            <field id="soiltemp_slope03"
     679                   long_name="Soil temperature for slope 03"
     680                   unit="K" />
     681            <field id="soiltemp_slope04"
     682                   long_name="Soil temperature for slope 04"
     683                   unit="K" />
     684            <field id="soiltemp_slope05"
     685                   long_name="Soil temperature for slope 05"
     686                   unit="K" />
     687            <field id="soiltemp_slope06"
     688                   long_name="Soil temperature for slope 06"
     689                   unit="K" />
     690            <field id="soiltemp_slope07"
     691                   long_name="Soil temperature for slope 07"
     692                   unit="K" />
    465693            <field id="inertiedat"
    466694                   long_name="Soil thermal inertia"
    467695                   unit="J/kg/K" />
     696            <field id="inertiesoil_slope01"
     697                   long_name="Soil thermal inertia for slope 01"
     698                   unit="J/kg/K" />
     699            <field id="inertiesoil_slope02"
     700                   long_name="Soil thermal inertia for slope 02"
     701                   unit="J/kg/K" />
     702            <field id="inertiesoil_slope03"
     703                   long_name="Soil thermal inertia for slope 03"
     704                   unit="J/kg/K" />
     705            <field id="inertiesoil_slope04"
     706                   long_name="Soil thermal inertia for slope 04"
     707                   unit="J/kg/K" />
     708            <field id="inertiesoil_slope05"
     709                   long_name="Soil thermal inertia for slope 05"
     710                   unit="J/kg/K" />
     711            <field id="inertiesoil_slope06"
     712                   long_name="Soil thermal inertia for slope 06"
     713                   unit="J/kg/K" />
     714            <field id="inertiesoil_slope07"
     715                   long_name="Soil thermal inertia for slope 07"
     716                   unit="J/kg/K" />
     717            <field id="Waterdensity_soil_slope01 for slope 01"
     718                   long_name="rhowater_soil"
     719                   unit="J/kg/K" />
     720            <field id="Waterdensity_soil_slope02 for slope 02"
     721                   long_name="rhowater_soil"
     722                   unit="J/kg/K" />
     723            <field id="Waterdensity_soil_slope03 for slope 03"
     724                   long_name="rhowater_soil"
     725                   unit="J/kg/K" />
     726            <field id="Waterdensity_soil_slope04 for slope 04"
     727                   long_name="rhowater_soil"
     728                   unit="J/kg/K" />
     729            <field id="Waterdensity_soil_slope05 for slope 05"
     730                   long_name="rhowater_soil"
     731                   unit="J/kg/K" />
     732            <field id="Waterdensity_soil_slope06 for slope 06"
     733                   long_name="rhowater_soil"
     734                   unit="J/kg/K" />
     735            <field id="Waterdensity_soil_slope07 for slope 07"
     736                   long_name="rhowater_soil"
     737                   unit="J/kg/K" />
     738
    468739        </field_group>
    469740
  • trunk/LMDZ.MARS/libf/phymars/call_dayperi.F

    r2448 r2963  
    1717c   Lsperi      solar longitude (Ls) of perohelion (rad)
    1818c   e_elips       Excentricity
     19c   real year_day      ! number of sols per Mars yar
    1920c
    2021c   output
    2122c   ------
    2223c   dayperi       Martian date at perihelion (sol)
    23 c   real year_day      ! number of sols per Mars yar
     24
    2425c-----------------------------------------------------------------------
    2526
  • trunk/LMDZ.MARS/libf/phymars/iniorbit.F

    r1382 r2963  
    6868
    6969      timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
     70      if(timeperi.lt.0) timeperi=-timeperi
    7071      PRINT*,'iniorbit: Perihelion solar long. Ls (deg)=',
    7172     &       360.-timeperi*180./pi
  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r2953 r2963  
    545545      logical :: write_restart
    546546
     547! Variable for ice table
     548      REAL :: rhowater_surf(ngrid,nslope)
     549      REAL :: rhowater_surf_sat(ngrid,nslope)
     550      REAL :: rhowater_soil(ngrid,nsoilmx,nslope)
     551      REAL,PARAMETER  :: alph_clap = -6143.7
     552      REAL,PARAMETER :: beta_clap = 28.9074
     553      REAL :: pvap_surf(ngrid)
     554      REAL,PARAMETER :: m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
     555      REAL,PARAMETER :: m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
     556      REAL :: ztmp1,ztmp2
     557
    547558c=======================================================================
    548559      pdq(:,:,:) = 0.
     
    789800      call update_xios_timestep     
    790801#endif
     802
     803
    791804
    792805
     
    14981511             dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:)
    14991512          ENDDO
     1513
    15001514          call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif,
    15011515     &   albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp)
     
    15591573           ENDDO
    15601574         ENDDO
     1575
    15611576         IF (turb_resolved) THEN
    15621577            write(*,*) 'Turbulent-resolving mode !'
     
    21602175      endif ! of if (callthermos)
    21612176#endif
     2177
    21622178c-----------------------------------------------------------------------
    21632179c   11. Carbon dioxide condensation-sublimation:
     
    38223838           endif !not.scavenging
    38233839        ENDIF ! of IF (water)
     3840
     3841! Output needed by the PEM
     3842          DO ig = 1,ngrid
     3843            ztmp1 =(1/m_co2 - 1/m_noco2)
     3844            ztmp2=1/m_noco2
     3845            pvap_surf(ig) = 1/(ztmp1*zq(ig,1,igcm_co2)+ztmp2)
     3846     &      * zq(ig,1,igcm_h2o_vap)/(mmol(igcm_h2o_vap)*1.e-3)*ps(ig)
     3847
     3848            DO islope = 1,nslope
     3849             rhowater_surf_sat(ig,islope)  =
     3850     &         exp(alph_clap/tsurf(ig,islope)+beta_clap)
     3851     &         / tsurf(ig,islope)
     3852 
     3853             if(qsurf(ig,igcm_h2o_ice,islope).gt.(1.e-4)) then
     3854               rhowater_surf(ig,islope) =
     3855     &         exp(alph_clap/tsurf(ig,islope)+beta_clap)
     3856     &         / tsurf(ig,islope)
     3857             else
     3858              rhowater_surf(ig,islope) = pvap_surf(ig)
     3859     &         / tsurf(ig,islope)
     3860             endif
     3861             DO isoil = 1,nsoilmx
     3862             rhowater_soil(ig,isoil,islope) =
     3863     &         exp(alph_clap/tsoil(ig,isoil,islope)+beta_clap)
     3864     &         / tsoil(ig,isoil,islope)
     3865             ENDDO
     3866          ENDDO
     3867        ENDDO
     3868
     3869      DO islope = 1,nslope
     3870        write(str2(1:2),'(i2.2)') islope
     3871        CALL send_xios_field("Waterdensity_soil_slope"//str2,
     3872     &     rhowater_soil(:,:,islope))
     3873        CALL send_xios_field("Waterdensity_surface"//str2,
     3874     &     rhowater_surf(:,islope))
     3875      ENDDO
     3876
     3877      CALL send_xios_field("h2o_layer1",zq(:,1,igcm_h2o_vap))
     3878      CALL send_xios_field("co2_layer1",zq(:,1,igcm_co2))
     3879
    38243880!PREVIOUSLY IN 1D ONLY
    38253881
  • trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90

    r2953 r2963  
    532532                          'lapse rate in the rocket dust storm', &
    533533                                    'K/m',lapserate(:,:))
    534         call write_output('rds_deltahr', &
    535                           'extra heating rate in the rocket dust storm', &
    536                                     'K/s',deltahr(:,:))
     534!        call write_output('rds_deltahr', &
     535!                          'extra heating rate in the rocket dust storm', &
     536!                                    'K/s',deltahr(:,:))
    537537!        call write_output('scheme','which scheme',&
    538538!                                                   ' ',real(scheme(:)))
  • trunk/LMDZ.MARS/util/startarchive2icosa/rearrange_startphy.f90

    r2941 r2963  
    270270    endif
    271271  enddo
     272
     273     print *, "All is well that ends well"
    272274 
    273275END PROGRAM rearrange_startphy
Note: See TracChangeset for help on using the changeset viewer.