Index: trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90 (revision 2962) +++ trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90 (revision 2963) @@ -25,39 +25,11 @@ !PEM parameters - - - year_bp_ini=0. - CALL getin('year_bp_ini', year_bp_ini) - - dt_pem=1 - CALL getin('dt_pem', dt_pem) - - water_ice_criterion=0.2 - CALL getin('water_ice_criterion', water_ice_criterion) - - co2_ice_criterion=0.2 - CALL getin('co2_ice_criterion', co2_ice_criterion) - - ps_criterion = 0.15 - CALL getin('ps_criterion',ps_criterion) +! #---------- ORBITAL parameters --------------# evol_orbit_pem=.false. CALL getin('evol_orbit_pem', evol_orbit_pem) - Max_iter_pem=99999999 - CALL getin('Max_iter_pem', Max_iter_pem) - - co2glaciersflow = .true. - CALL getin('co2glaciersflow', co2glaciersflow) - - soil_pem=.true. - CALL getin('soil_pem', soil_pem) - - adsorption_pem = .true. - CALL getin('adsorption_pem',adsorption_pem) - - fluxgeo = 0. - CALL getin('Fluxgeo_PEM',fluxgeo) - print*,'Flux Geothermal is set to',fluxgeo + year_bp_ini=0. + CALL getin('year_bp_ini', year_bp_ini) var_obl = .true. @@ -72,4 +44,42 @@ CALL getin('var_lsp',var_lsp) print*,'Does Ls peri vary ?',var_lsp + +! #---------- Stopping criterion parameters --------------# + + Max_iter_pem=99999999 + CALL getin('Max_iter_pem', Max_iter_pem) + + water_ice_criterion=0.2 + CALL getin('water_ice_criterion', water_ice_criterion) + + co2_ice_criterion=0.2 + CALL getin('co2_ice_criterion', co2_ice_criterion) + + ps_criterion = 0.15 + CALL getin('ps_criterion',ps_criterion) + + dt_pem=1 + CALL getin('dt_pem', dt_pem) + +! #---------- Subsurface parameters --------------# + + soil_pem=.true. + CALL getin('soil_pem', soil_pem) + + adsorption_pem = .true. + CALL getin('adsorption_pem',adsorption_pem) + + co2glaciersflow = .true. + CALL getin('co2glaciersflow', co2glaciersflow) + + reg_thprop_dependp = .false. + CALL getin('reg_thprop_dependp',reg_thprop_dependp) + print*, 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp + +! #---------- Layering parameters --------------# + + fluxgeo = 0. + CALL getin('Fluxgeo_PEM',fluxgeo) + print*,'Flux Geothermal is set to',fluxgeo depth_breccia = 10. @@ -80,8 +90,4 @@ CALL getin('depth_bedrock',depth_bedrock) print*,'Depth of bedrock is set to',depth_bedrock - - reg_thprop_dependp = .false. - CALL getin('reg_thprop_dependp',reg_thprop_dependp) - print*, 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp icetable_equilibrium = .true. Index: trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90 (revision 2962) +++ trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90 (revision 2963) @@ -86,5 +86,4 @@ enddo endif - endif negative_part = 0. Index: trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90 (revision 2962) +++ trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90 (revision 2963) @@ -46,4 +46,5 @@ logical obl_not_found, ex_not_found,lsp_not_found !Loop variable (first call) + logical max_lsp_modulo,min_lsp_modulo ! ********************************************************************** @@ -109,4 +110,9 @@ min_lsp=timeperi_ls-max_change_lsp + max_lsp_modulo=.false. + min_lsp_modulo=.false. + if(max_lsp.gt.360) max_lsp_modulo=.true. + if(min_lsp.lt.0) min_lsp_modulo=.true. + !End Constant max change case @@ -133,30 +139,32 @@ max_lsp_iter=999999999999 - do ilask=last_ilask+1,1,-1 + do ilask=last_ilask,1,-1 if((oblask(ilask).GT.max_obl).and. obl_not_found ) then - max_obl_iter=(max_obl-oblask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & - / (oblask(ilask-1)-oblask(ilask)) + max_obl_iter=(max_obl-oblask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) & + / (oblask(ilask)-oblask(ilask-1))+yearlask(ilask-1)-Year obl_not_found=.FALSE. elseif((oblask(ilask).LT.min_obl).and. obl_not_found ) then - max_obl_iter=(min_obl-oblask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & - / (oblask(ilask-1)-oblask(ilask)) + max_obl_iter=(min_obl-oblask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) & + / (oblask(ilask)-oblask(ilask-1))+yearlask(ilask-1)-Year obl_not_found=.FALSE. endif if((exlask(ilask).GT.max_ex).and. ex_not_found ) then - max_ex_iter=(max_ex-exlask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & - / (exlask(ilask-1)-exlask(ilask)) + max_ex_iter=(max_ex-exlask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) & + / (exlask(ilask)-exlask(ilask-1))+yearlask(ilask-1)-Year ex_not_found=.FALSE. elseif((exlask(ilask).LT.min_ex ).and. ex_not_found ) then - max_ex_iter=(min_ex-exlask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & - / (exlask(ilask-1)-exlask(ilask)) + max_ex_iter=(min_ex-exlask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) & + / (exlask(ilask)-exlask(ilask-1))+yearlask(ilask-1)-Year ex_not_found=.FALSE. endif + if(max_lsp_modulo .and. lsplask(ilask).lt.180) lsplask(ilask)=lsplask(ilask)+360. + if(min_lsp_modulo .and. lsplask(ilask).gt.180) lsplask(ilask)=lsplask(ilask)-360. if((lsplask(ilask).GT.max_lsp).and. lsp_not_found ) then - max_lsp_iter=(max_lsp-lsplask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & - / (lsplask(ilask-1)-lsplask(ilask)) + max_lsp_iter=(max_lsp-lsplask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) & + / (lsplask(ilask)-lsplask(ilask-1))+yearlask(ilask-1)-Year lsp_not_found=.FALSE. elseif((lsplask(ilask).LT.min_lsp ).and. lsp_not_found ) then - max_lsp_iter=(min_lsp-lsplask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & - / (lsplask(ilask-1)-lsplask(ilask)) + max_lsp_iter=(min_lsp-lsplask(ilask-1)) * (yearlask(ilask)-yearlask(ilask-1)) & + / (lsplask(ilask)-lsplask(ilask-1))+yearlask(ilask-1)-Year lsp_not_found=.FALSE. endif @@ -171,6 +179,6 @@ print *, "Maximum number of iteration for the ex. parameter=", max_ex_iter - print *, "Maximum lsp accepted=", max_obl - print *, "Minimum lsp accepted=", min_obl + print *, "Maximum lsp accepted=", max_lsp + print *, "Minimum lsp accepted=", min_lsp print *, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 2962) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 2963) @@ -32,6 +32,6 @@ !module needed for INITIALISATION #ifndef CPP_STD - use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa - use surfdat_h, only: tsurf, co2ice, emis,& + use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa,inertiesoil + use surfdat_h, only: tsurf, emis,& qsurf,watercap, ini_surfdat_h, & albedodat, zmea, zstd, zsig, zgam, zthe, & @@ -40,10 +40,10 @@ use dimradmars_mod, only: totcloudfrac, albedo use dust_param_mod, only: tauscaling - use tracer_mod, only: noms,igcm_h2o_ice ! tracer names + use tracer_mod, only: noms,igcm_h2o_ice,igcm_co2 ! tracer names #else use comsoil_h, only: nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe, & emissiv - use tracer_h, only: noms,igcm_h2o_ice,igcm_co2_ice ! tracer names + use tracer_h, only: noms,igcm_h2o_ice,igcm_co2 ! tracer names use phys_state_var_mod #endif @@ -73,10 +73,6 @@ USE mod_const_mpi, ONLY: COMM_LMDZ USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, & - subslope_dist,co2ice_slope, & - tsurf_slope,tsoil_slope,fluxgrd_slope,& - fluxrad_sky_slope,sky_slope,callsubslope,& - co2iceflow, beta_slope, capcal_slope,& - albedo_slope,emiss_slope,qsurf_slope,& - iflat,major_slope,ini_comslope_h,fluxgeo_slope + subslope_dist,iflat, & + major_slope,ini_comslope_h use time_phylmdz_mod, only: daysec,dtphys USE comconst_mod, ONLY: rad,g,r,cpp,pi @@ -94,9 +90,9 @@ TI_PEM,inertiedat_PEM, & ! soil thermal inertia tsoil_PEM, mlayer_PEM,layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths - fluxgeo, & ! geothermal flux for the PEM and GCM + fluxgeo, & ! geothermal flux for the PEM and GCM water_reservoir ! Water ressources - use adsorption_mod, only : regolith_adsorption,adsorption_pem, & ! bool to check if adsorption, main subroutine - ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! allocate arrays - co2_adsorbded_phys, h2o_adsorbded_phys ! mass of co2 and h2O adsorbded + use adsorption_mod, only : regolith_adsorption,adsorption_pem, & ! bool to check if adsorption, main subroutine + ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! allocate arrays + co2_adsorbded_phys, h2o_adsorbded_phys ! mass of co2 and h2O adsorbded !!! For orbit parameters @@ -198,14 +194,13 @@ !!!!!!!!!!!!!!!!!!!!!!!! SLOPE - REAL ,allocatable :: watercap_slope(:,:) ! Physics x Nslope: watercap per slope - REAL :: watercap_slope_saved ! Value saved at the previous time step + REAL :: watercap_saved ! Value saved at the previous time step REAL , dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2] REAL , dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2] REAL , dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field : minimum of water ice at each point for the first year [kg/m^2] REAL , dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field : minimum of water ice at each point for the second year [kg/m^2] - REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_slope ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] - REAL, dimension(:,:),allocatable :: initial_co2_ice_sublim_slope ! physical point field : Logical array indicating sublimating point of co2 ice - REAL, dimension(:,:),allocatable :: initial_h2o_ice_slope ! physical point field : Logical array indicating if there is water ice at initial state - REAL, dimension(:,:),allocatable :: initial_co2_ice_slope ! physical point field : Logical array indicating if there is co2 ice at initial state + REAL , dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] + REAL, dimension(:,:),allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice + REAL, dimension(:,:),allocatable :: initial_h2o_ice ! physical point field : Logical array indicating if there is water ice at initial state + REAL, dimension(:,:),allocatable :: initial_co2_ice ! physical point field : Logical array indicating if there is co2 ice at initial state REAL, dimension(:,:),allocatable :: tendencies_co2_ice ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year 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 @@ -225,6 +220,5 @@ REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] - REAL, ALLOCATABLE :: inertiesoil(:,:) !Physic x Depth Thermal inertia of the mesh for restart [SI] - REAL, ALLOCATABLE :: TI_GCM(:,:,:) ! Same but for the start + REAL,ALLOCATABLE :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] REAL,ALLOCATABLE :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K] @@ -310,28 +304,4 @@ !----------------------------Initialisation : READ some constant of startfi_evol.nc --------------------- -! In the gcm, these values are given to the physic by the dynamic. -! Here we simply read them in the startfi_evol.nc file - status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) - - allocate(longitude(ngrid)) - allocate(latitude(ngrid)) - allocate(cell_area(ngrid)) - - status = nf90_inq_varid(ncid, "longitude", lonvarid) - status = nf90_get_var(ncid, lonvarid, longitude) - - status = nf90_inq_varid(ncid, "latitude", latvarid) - status = nf90_get_var(ncid, latvarid, latitude) - - status = nf90_inq_varid(ncid, "area", areavarid) - status = nf90_get_var(ncid, areavarid, cell_area) - - call ini_comsoil_h(ngrid) - - status = nf90_inq_varid(ncid, "soildepth", sdvarid) - status = nf90_get_var(ncid, sdvarid, mlayer) - - status =nf90_close(ncid) - !----------------------------READ start.nc --------------------- @@ -360,11 +330,29 @@ iflag_phys) +! In the gcm, these values are given to the physic by the dynamic. +! Here we simply read them in the startfi_evol.nc file + status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) + + allocate(longitude(ngrid)) + allocate(latitude(ngrid)) + allocate(cell_area(ngrid)) + + status = nf90_inq_varid(ncid, "longitude", lonvarid) + status = nf90_get_var(ncid, lonvarid, longitude) + + status = nf90_inq_varid(ncid, "latitude", latvarid) + status = nf90_get_var(ncid, latvarid, latitude) + + status = nf90_inq_varid(ncid, "area", areavarid) + status = nf90_get_var(ncid, areavarid, cell_area) + + status = nf90_inq_varid(ncid, "soildepth", sdvarid) + status = nf90_get_var(ncid, sdvarid, mlayer) + + status =nf90_close(ncid) + !----------------------------READ startfi.nc --------------------- ! First we read the initial state (starfi.nc) - - allocate(watercap_slope(ngrid,nslope)) - allocate(TI_GCM(ngrid,nsoilmx,nslope)) - allocate(inertiesoil(ngrid,nsoilmx)) #ifndef CPP_STD @@ -373,11 +361,6 @@ day_ini,time_phys, & tsurf,tsoil,albedo,emis, & - q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar, & - watercap,inertiesoil,nslope,tsurf_slope, & - tsoil_slope,co2ice_slope,def_slope,def_slope_mean, & - subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM, & - qsurf_slope,watercap_slope) - - + q2,qsurf,tauscaling,totcloudfrac,wstar, & + watercap,def_slope,def_slope_mean,subslope_dist) ! Remove unphysical values of surface tracer @@ -385,6 +368,6 @@ DO nnq=1,nqtot DO islope=1,nslope - if(qsurf_slope(i,nnq,islope).LT.0) then - qsurf_slope(i,nnq,islope)=0. + if(qsurf(i,nnq,islope).LT.0) then + qsurf(i,nnq,islope)=0. endif enddo @@ -409,5 +392,5 @@ allocate(co2ice(ngrid)) - co2ice(:)=qsurf(:,igcm_co2_ice) + co2ice(:)=qsurf(:,igcm_co2) co2ice_slope(:,1)=co2ice(:) tsurf_slope(:,1)=tsurf(:) @@ -433,4 +416,5 @@ DO nnq=1,nqtot if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq + if(noms(nnq).eq."co2") igcm_co2 = nnq ENDDO @@ -492,5 +476,5 @@ allocate(q_co2_PEM_phys(ngrid,timelen)) allocate(q_h2o_PEM_phys(ngrid,timelen)) - allocate(co2_ice_GCM_slope(ngrid,nslope,timelen)) + allocate(co2_ice_GCM(ngrid,nslope,timelen)) allocate(watersurf_density_ave(ngrid,nslope)) allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) @@ -508,5 +492,5 @@ 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,& - tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM_slope, & + tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, & watersurf_density_ave,watersoil_density_timeseries) @@ -517,5 +501,5 @@ 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, & - tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM_slope, & + tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, & watersurf_density_ave,watersoil_density_timeseries) @@ -543,5 +527,5 @@ if(soil_pem) then - call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM,TI_PEM) + call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) DO l=1,nsoilmx tsoil_ave_PEM_yr1(:,l,:)=tsoil_ave_yr1(:,l,:) @@ -605,7 +589,7 @@ !---------------------------- Save initial PCM situation --------------------- - allocate(initial_co2_ice_sublim_slope(ngrid,nslope)) - allocate(initial_co2_ice_slope(ngrid,nslope)) - allocate(initial_h2o_ice_slope(ngrid,nslope)) + allocate(initial_co2_ice_sublim(ngrid,nslope)) + allocate(initial_co2_ice(ngrid,nslope)) + allocate(initial_h2o_ice(ngrid,nslope)) ! We save the places where water ice is sublimating @@ -618,19 +602,19 @@ do islope=1,nslope if (tendencies_co2_ice(i,islope).LT.0) then - initial_co2_ice_sublim_slope(i,islope)=1. + initial_co2_ice_sublim(i,islope)=1. ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope) else - initial_co2_ice_sublim_slope(i,islope)=0. + initial_co2_ice_sublim(i,islope)=0. endif - if (co2ice_slope(i,islope).GT.0) then - initial_co2_ice_slope(i,islope)=1. + if (qsurf(i,igcm_co2,islope).GT.0) then + initial_co2_ice(i,islope)=1. else - initial_co2_ice_slope(i,islope)=0. + initial_co2_ice(i,islope)=0. endif if (tendencies_h2o_ice(i,islope).LT.0) then - initial_h2o_ice_slope(i,islope)=1. + initial_h2o_ice(i,islope)=1. ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope) else - initial_h2o_ice_slope(i,islope)=0. + initial_h2o_ice(i,islope)=0. endif enddo @@ -674,5 +658,5 @@ call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_PEM_yr1,tsoil_PEM,porefillingice_depth,porefillingice_thickness, & tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,& - tendencies_h2o_ice,tendencies_co2_ice,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,& + tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,& watersurf_density_ave,watersoil_density_PEM_ave, & co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir) @@ -680,5 +664,5 @@ do ig = 1,ngrid do islope = 1,nslope - 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.) + 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.) enddo enddo @@ -829,14 +813,14 @@ print *, "Evolution of h2o ice" - call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope) + call evol_h2o_ice_s_slope(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope) DO islope=1, nslope write(str2(1:2),'(i2.2)') islope - call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf_slope(:,igcm_h2o_ice,islope)) + call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope)) call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope)) ENDDO print *, "Evolution of co2 ice" - call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope) + call evol_co2_ice_s_slope(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope) !------------------------ @@ -853,5 +837,5 @@ if (co2glaciersflow) then call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,& - global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh) + global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh) endif @@ -860,5 +844,5 @@ DO islope=1, nslope write(str2(1:2),'(i2.2)') islope - call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,co2ice_slope(:,islope)) + call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) ENDDO @@ -881,9 +865,9 @@ DO ig = 1,ngrid DO islope = 1,nslope - 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 + 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 if(latitude_deg(ig).gt.0) then do ig_loop=ig,ngrid DO islope_loop=islope,iflat,-1 - if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then + if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop) bool_sublim=1 @@ -898,5 +882,5 @@ do ig_loop=ig,1,-1 DO islope_loop=islope,iflat - if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then + if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop) bool_sublim=1 @@ -909,17 +893,17 @@ enddo endif - initial_co2_ice_slope(ig,islope)=0 - if ((co2ice_slope(ig,islope).lt.1e-10).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then - albedo_slope(ig,1,islope) = albedo_h2o_frost - albedo_slope(ig,2,islope) = albedo_h2o_frost + initial_co2_ice(ig,islope)=0 + if ((qsurf(ig,igcm_co2,islope).lt.1e-10).and. (qsurf(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then + albedo(ig,1,islope) = albedo_h2o_frost + albedo(ig,2,islope) = albedo_h2o_frost else - albedo_slope(ig,1,islope) = albedodat(ig) - albedo_slope(ig,2,islope) = albedodat(ig) - emiss_slope(ig,islope) = emissiv + albedo(ig,1,islope) = albedodat(ig) + albedo(ig,2,islope) = albedodat(ig) + emis(ig,islope) = emissiv endif - elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2 + elseif( (qsurf(ig,igcm_co2,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2 ave=0. do t=1,timelen - if(co2_ice_GCM_slope(ig,islope,t).gt.1e-3) then + if(co2_ice_GCM(ig,islope,t).gt.1e-3) then ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.)) else @@ -938,5 +922,5 @@ do ig = 1,ngrid do islope = 1,nslope - tsurf_slope(ig,islope) = tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope)) + tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope)) enddo enddo @@ -944,5 +928,5 @@ DO islope=1, nslope write(str2(1:2),'(i2.2)') islope - call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_slope(:,islope)) + call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) ENDDO @@ -990,10 +974,10 @@ ! II_d.4 Update the soil thermal properties - call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, & + call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, & porefillingice_depth,porefillingice_thickness,TI_PEM) ! II_d.5 Update the mass of the regolith adsorbded if(adsorption_pem) then - call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf_slope(:,igcm_h2o_ice,:),co2ice_slope, & + call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:), & tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded) @@ -1013,5 +997,5 @@ print *, "Adaptation of the new co2 tendencies to the current pressure" - call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,co2ice_slope,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,& + call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,& global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope) @@ -1027,8 +1011,8 @@ !------------------------ - call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope) - - call criterion_co2_stop(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,STOPPING_pressure,ngrid, & - initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope) + call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice) + + call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, & + initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope) year_iter=year_iter+dt_pem @@ -1084,4 +1068,5 @@ ! III_a.1 Ice update (for startfi) +#ifdef CPP_STD ! Co2 ice DO ig = 1,ngrid @@ -1092,8 +1077,7 @@ cos(pi*def_slope_mean(islope)/180.) ENDDO -#ifdef CPP_STD - qsurf(ig,igcm_co2_ice)=co2ice(ig) + qsurf(ig,igcm_co2)=co2ice(ig) + ENDDO ! of DO ig=1,ngrid #endif - ENDDO ! of DO ig=1,ngrid ! H2O ice @@ -1102,13 +1086,13 @@ watercap_sum=0. DO islope=1,nslope - watercap_slope_saved = watercap_slope(ig,islope) - if(qsurf_slope(ig,igcm_h2o_ice,islope).GT. (watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then - 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.)) + watercap_saved = watercap(ig,islope) + if(qsurf(ig,igcm_h2o_ice,islope).GT. (watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then + 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.)) else - 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.)) - qsurf_slope(ig,igcm_h2o_ice,islope)=0. + 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.)) + qsurf(ig,igcm_h2o_ice,islope)=0. endif - watercap_sum=watercap_sum+(watercap_slope(ig,islope)-watercap_slope_saved)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) - watercap_slope(ig,islope)=0. + watercap_sum=watercap_sum+(watercap(ig,islope)-watercap_saved)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) + watercap(ig,islope)=0. enddo water_reservoir(ig)=water_reservoir(ig)+watercap_sum @@ -1116,18 +1100,9 @@ enddo - DO ig = 1,ngrid - qsurf(ig,igcm_h2o_ice) = 0. - DO islope = 1,nslope - qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) & - * subslope_dist(ig,islope) / & - cos(pi*def_slope_mean(islope)/180.) - ENDDO - ENDDO ! of DO ig=1,ngrid - DO ig=1,ngrid DO islope=1,nslope - if(qsurf_slope(ig,igcm_h2o_ice,islope).GT.500) then + if(qsurf(ig,igcm_h2o_ice,islope).GT.500) then watercaptag(ig)=.true. - qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-250 + qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)-250 water_reservoir(ig)=water_reservoir(ig)+250 endif @@ -1136,74 +1111,18 @@ DO ig=1,ngrid - if(water_reservoir(ig).LT. 10) then - watercaptag(ig)=.false. - qsurf(ig,igcm_h2o_ice)=qsurf(ig,igcm_h2o_ice)+water_reservoir(ig) - DO islope=1,nslope - qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+water_reservoir(ig) - ENDDO - endif - enddo - - DO ig = 1,ngrid - watercap(ig) = 0. - DO islope = 1,nslope - watercap(ig) = watercap(ig) + watercap_slope(ig,islope) & - * subslope_dist(ig,islope) / & - cos(pi*def_slope_mean(islope)/180.) - ENDDO - ENDDO ! of DO ig=1,ngrid + if(water_reservoir(ig).LT. 10) then + watercaptag(ig)=.false. + DO islope=1,nslope + qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+water_reservoir(ig) + ENDDO + endif + enddo ! III_a.2 Tsoil update (for startfi) if(soil_pem) then - fluxgeo_slope(:,:) = fluxgeo - call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM) - tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen) + call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) + tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen) endif !soil_pem - -#ifndef CPP_STD - DO ig = 1,ngrid - DO iloop = 1,nsoilmx - tsoil(ig,iloop) = 0. - inertiesoil(ig,iloop) = 0. - DO islope = 1,nslope - tsoil(ig,iloop) = tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) & - * subslope_dist(ig,islope) - inertiesoil(ig,iloop) = inertiesoil(ig,iloop) + TI_GCM(ig,iloop,islope) & - * subslope_dist(ig,islope) - ENDDO - ENDDO - ENDDO ! of DO ig=1,ngrid - -! III_a.3 Surface optical properties (for startfi) - - DO ig = 1,ngrid - DO l = 1,2 - albedo(ig,l) =0. - DO islope = 1,nslope - albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) & - *subslope_dist(ig,islope) - ENDDO - ENDDO - ENDDO - - DO ig = 1,ngrid - emis(ig) =0. - DO islope = 1,nslope - emis(ig)= emis(ig)+emiss_slope(ig,islope) & - *subslope_dist(ig,islope) - ENDDO - ENDDO - - DO ig = 1,ngrid - tsurf(ig) = 0. - DO islope = 1,nslope - tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)*tsurf_slope(ig,islope))**4 & - * subslope_dist(ig,islope) - ENDDO - tsurf(ig) = tsurf(ig)**(0.25)/emis(ig) - ENDDO ! of DO ig=1,ngrid - -#endif ! III_a.4 Pressure (for start) @@ -1250,5 +1169,5 @@ DO ig=1,ngrid DO l=1,llm-1 - if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") ) then + 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 extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1)) q(ig,l,nnq)=1. @@ -1295,15 +1214,12 @@ nsoilmx,ngrid,nlayer,nq, & ptimestep,pday,0.,cell_area, & - albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, & - hmons,summit,base,nslope,def_slope, & + albedodat,inertiedat,def_slope, & subslope_dist) call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & ptimestep,ztime_fin, & - tsurf,tsoil,co2ice,albedo,emis, & - q2,qsurf,tauscaling,totcloudfrac,wstar, & - watercap,inertiesoil,nslope,co2ice_slope, & - tsurf_slope,tsoil_slope, albedo_slope, & - emiss_slope,qsurf_slope,watercap_slope, TI_GCM) + tsurf,tsoil,inertiesoil,albedo, & + emis,q2,qsurf,tauscaling,totcloudfrac,wstar,& + watercap) #else call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & @@ -1346,5 +1262,5 @@ deallocate(q_co2_PEM_phys) deallocate(q_h2o_PEM_phys) - deallocate(co2_ice_GCM_slope) + deallocate(co2_ice_GCM) deallocate(watersurf_density_ave) deallocate(watersoil_density_timeseries) Index: trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 (revision 2962) +++ trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 (revision 2963) @@ -91,10 +91,10 @@ print *, "Downloading data for vmr co2..." - CALL get_var3("co2_cropped" ,q_co2_dyn) + CALL get_var3("co2_layer1" ,q_co2_dyn) print *, "Downloading data for vmr co2 done" print *, "Downloading data for vmr h20..." - CALL get_var3("h2o_cropped" ,q_h2o_dyn) + CALL get_var3("h2o_layer1" ,q_h2o_dyn) print *, "Downloading data for vmr h2o done" @@ -143,12 +143,12 @@ if(soil_pem) then - print *, "Downloading data for tsoil_slope ..." - -DO islope=1,nslope - write(num,fmt='(i2.2)') islope - call get_var4("tsoil_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:)) -ENDDO - - print *, "Downloading data for tsoil_slope done" + print *, "Downloading data for soiltemp_slope ..." + +DO islope=1,nslope + write(num,fmt='(i2.2)') islope + call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:)) +ENDDO + + print *, "Downloading data for soiltemp_slope done" print *, "Downloading data for watersoil_density ..." @@ -182,5 +182,5 @@ if(soil_pem) then - call get_var4("tsoil",tsoil_gcm_dyn(:,:,:,1,:)) + call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:)) endif !soil_pem endif !nslope=1 Index: trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 (revision 2962) +++ trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 (revision 2963) @@ -16,5 +16,5 @@ #ifndef CPP_STD USE comconst_mod, ONLY: pi - USE planete_h, ONLY: e_elips, obliquit, timeperi + USE planete_h, ONLY: e_elips, obliquit, timeperi, periheli,aphelie,p_elips,peri_day,year_day #else use planete_mod, only: e_elips, obliquit, timeperi @@ -45,4 +45,6 @@ real :: timeperi_ls ! time of peri in ls integer ilask ! Loop variable + real :: halfaxe ! Million of km + real :: unitastr @@ -73,10 +75,24 @@ endif if(var_lsp) then + if(lsplask(ilask)-lsplask(ilask+1) .gt.200) then + if(lsplask(ilask).lt.lsplask(ilask+1)) then + lsplask(ilask)=lsplask(ilask)+360 + else + lsplask(ilask+1)=lsplask(ilask+1)+360 + endif + endif timeperi_ls=lsplask(ilask+1)+(lsplask(ilask)-lsplask(ilask+1))*(Year-yearlask(ilask+1))/(yearlask(ilask)-yearlask(ilask+1)) + timeperi_ls=modulo(timeperi_ls,360.) endif exit endif enddo + halfaxe=227.94 timeperi=timeperi_ls*2*pi/360 + periheli = halfaxe*(1-e_elips) + aphelie = halfaxe*(1+e_elips) + unitastr=149.597927 + p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr + call call_dayperi(timeperi,e_elips,peri_day,year_day) print *, "recomp_orb_param, Final year of the PEM run:", Year Index: trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90 (revision 2962) +++ trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90 (revision 2963) @@ -51,8 +51,4 @@ integer :: numyear - - - - DO numyear=1, 2 write(*,*) 'numyear',numyear @@ -61,8 +57,8 @@ !integer :: ncid, status !... -status = nf90_open(path = "Xdiurnalave"//str2//".nc", mode = nf90_nowrite, ncid = ncid1) -if(status /= nf90_noerr) call handle_err(status) - -status = nf90_create(path = "Xdiurnalave"//str2//".nc_new", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2) +status = nf90_open(path = "data2reshape"//str2//".nc", mode = nf90_nowrite, ncid = ncid1) +if(status /= nf90_noerr) call handle_err(status) + +status = nf90_create(path = "datareshaped"//str2//".nc", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2) if(status /= nf90_noerr) call handle_err(status)