Ignore:
Timestamp:
May 11, 2023, 5:40:03 PM (2 years 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/LMDZ.COMMON/libf/evolution
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified 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.

  • TabularUnified 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

  • TabularUnified 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

  • TabularUnified 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)

  • TabularUnified 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

  • TabularUnified 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

  • TabularUnified 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
Note: See TracChangeset for help on using the changeset viewer.