Changeset 2963
- Timestamp:
- May 11, 2023, 5:40:03 PM (20 months ago)
- Location:
- trunk
- Files:
-
- 5 added
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90
r2961 r2963 25 25 !PEM parameters 26 26 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 --------------# 43 28 44 29 evol_orbit_pem=.false. 45 30 CALL getin('evol_orbit_pem', evol_orbit_pem) 46 31 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) 62 34 63 35 var_obl = .true. … … 72 44 CALL getin('var_lsp',var_lsp) 73 45 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 74 84 75 85 depth_breccia = 10. … … 80 90 CALL getin('depth_bedrock',depth_bedrock) 81 91 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_dependp86 92 87 93 icetable_equilibrium = .true. -
trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90
r2950 r2963 86 86 enddo 87 87 endif 88 endif89 88 negative_part = 0. 90 89 -
trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90
r2894 r2963 46 46 47 47 logical obl_not_found, ex_not_found,lsp_not_found !Loop variable (first call) 48 logical max_lsp_modulo,min_lsp_modulo 48 49 49 50 ! ********************************************************************** … … 109 110 min_lsp=timeperi_ls-max_change_lsp 110 111 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 111 117 !End Constant max change case 112 118 … … 133 139 max_lsp_iter=999999999999 134 140 135 do ilask=last_ilask +1,1,-1141 do ilask=last_ilask,1,-1 136 142 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 139 145 obl_not_found=.FALSE. 140 146 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 143 149 obl_not_found=.FALSE. 144 150 endif 145 151 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 148 154 ex_not_found=.FALSE. 149 155 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 152 158 ex_not_found=.FALSE. 153 159 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. 154 162 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 157 165 lsp_not_found=.FALSE. 158 166 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 161 169 lsp_not_found=.FALSE. 162 170 endif … … 171 179 print *, "Maximum number of iteration for the ex. parameter=", max_ex_iter 172 180 173 print *, "Maximum lsp accepted=", max_ obl174 print *, "Minimum lsp accepted=", min_ obl181 print *, "Maximum lsp accepted=", max_lsp 182 print *, "Minimum lsp accepted=", min_lsp 175 183 print *, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter 176 184 -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r2961 r2963 32 32 !module needed for INITIALISATION 33 33 #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,& 36 36 qsurf,watercap, ini_surfdat_h, & 37 37 albedodat, zmea, zstd, zsig, zgam, zthe, & … … 40 40 use dimradmars_mod, only: totcloudfrac, albedo 41 41 use dust_param_mod, only: tauscaling 42 use tracer_mod, only: noms,igcm_h2o_ice ! tracer names42 use tracer_mod, only: noms,igcm_h2o_ice,igcm_co2 ! tracer names 43 43 #else 44 44 use comsoil_h, only: nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa 45 45 use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe, & 46 46 emissiv 47 use tracer_h, only: noms,igcm_h2o_ice,igcm_co2 _ice! tracer names47 use tracer_h, only: noms,igcm_h2o_ice,igcm_co2 ! tracer names 48 48 use phys_state_var_mod 49 49 #endif … … 73 73 USE mod_const_mpi, ONLY: COMM_LMDZ 74 74 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 81 77 use time_phylmdz_mod, only: daysec,dtphys 82 78 USE comconst_mod, ONLY: rad,g,r,cpp,pi … … 94 90 TI_PEM,inertiedat_PEM, & ! soil thermal inertia 95 91 tsoil_PEM, mlayer_PEM,layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths 96 fluxgeo, & ! geothermal flux for the PEM and GCM92 fluxgeo, & ! geothermal flux for the PEM and GCM 97 93 water_reservoir ! Water ressources 98 use adsorption_mod, only : regolith_adsorption,adsorption_pem, & ! bool to check if adsorption, main subroutine99 ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! allocate arrays100 co2_adsorbded_phys, h2o_adsorbded_phys ! mass of co2 and h2O adsorbded94 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 101 97 102 98 !!! For orbit parameters … … 198 194 199 195 !!!!!!!!!!!!!!!!!!!!!!!! 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 202 197 REAL , dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2] 203 198 REAL , dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2] 204 199 REAL , dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field : minimum of water ice at each point for the first year [kg/m^2] 205 200 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 ice208 REAL, dimension(:,:),allocatable :: initial_h2o_ice _slope! physical point field : Logical array indicating if there is water ice at initial state209 REAL, dimension(:,:),allocatable :: initial_co2_ice _slope! physical point field : Logical array indicating if there is co2 ice at initial state201 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 210 205 REAL, dimension(:,:),allocatable :: tendencies_co2_ice ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year 211 206 REAL, dimension(:,:),allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM … … 225 220 REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] 226 221 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 229 223 REAL,ALLOCATABLE :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] 230 224 REAL,ALLOCATABLE :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K] … … 310 304 !----------------------------Initialisation : READ some constant of startfi_evol.nc --------------------- 311 305 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 file314 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 336 306 !----------------------------READ start.nc --------------------- 337 307 … … 360 330 iflag_phys) 361 331 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 362 354 !----------------------------READ startfi.nc --------------------- 363 355 364 356 ! 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))369 357 370 358 #ifndef CPP_STD … … 373 361 day_ini,time_phys, & 374 362 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) 382 365 383 366 ! Remove unphysical values of surface tracer … … 385 368 DO nnq=1,nqtot 386 369 DO islope=1,nslope 387 if(qsurf _slope(i,nnq,islope).LT.0) then388 qsurf _slope(i,nnq,islope)=0.370 if(qsurf(i,nnq,islope).LT.0) then 371 qsurf(i,nnq,islope)=0. 389 372 endif 390 373 enddo … … 409 392 410 393 allocate(co2ice(ngrid)) 411 co2ice(:)=qsurf(:,igcm_co2 _ice)394 co2ice(:)=qsurf(:,igcm_co2) 412 395 co2ice_slope(:,1)=co2ice(:) 413 396 tsurf_slope(:,1)=tsurf(:) … … 433 416 DO nnq=1,nqtot 434 417 if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq 418 if(noms(nnq).eq."co2") igcm_co2 = nnq 435 419 ENDDO 436 420 … … 492 476 allocate(q_co2_PEM_phys(ngrid,timelen)) 493 477 allocate(q_h2o_PEM_phys(ngrid,timelen)) 494 allocate(co2_ice_GCM _slope(ngrid,nslope,timelen))478 allocate(co2_ice_GCM(ngrid,nslope,timelen)) 495 479 allocate(watersurf_density_ave(ngrid,nslope)) 496 480 allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) … … 508 492 509 493 call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1,& 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, & 511 495 watersurf_density_ave,watersoil_density_timeseries) 512 496 … … 517 501 518 502 call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, & 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, & 520 504 watersurf_density_ave,watersoil_density_timeseries) 521 505 … … 543 527 544 528 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) 546 530 DO l=1,nsoilmx 547 531 tsoil_ave_PEM_yr1(:,l,:)=tsoil_ave_yr1(:,l,:) … … 605 589 !---------------------------- Save initial PCM situation --------------------- 606 590 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)) 610 594 611 595 ! We save the places where water ice is sublimating … … 618 602 do islope=1,nslope 619 603 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. 621 605 ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope) 622 606 else 623 initial_co2_ice_sublim _slope(i,islope)=0.607 initial_co2_ice_sublim(i,islope)=0. 624 608 endif 625 if ( co2ice_slope(i,islope).GT.0) then626 initial_co2_ice _slope(i,islope)=1.609 if (qsurf(i,igcm_co2,islope).GT.0) then 610 initial_co2_ice(i,islope)=1. 627 611 else 628 initial_co2_ice _slope(i,islope)=0.612 initial_co2_ice(i,islope)=0. 629 613 endif 630 614 if (tendencies_h2o_ice(i,islope).LT.0) then 631 initial_h2o_ice _slope(i,islope)=1.615 initial_h2o_ice(i,islope)=1. 632 616 ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope) 633 617 else 634 initial_h2o_ice _slope(i,islope)=0.618 initial_h2o_ice(i,islope)=0. 635 619 endif 636 620 enddo … … 674 658 call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_PEM_yr1,tsoil_PEM,porefillingice_depth,porefillingice_thickness, & 675 659 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,& 677 661 watersurf_density_ave,watersoil_density_PEM_ave, & 678 662 co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir) … … 680 664 do ig = 1,ngrid 681 665 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.) 683 667 enddo 684 668 enddo … … 829 813 print *, "Evolution of h2o ice" 830 814 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) 832 816 833 817 DO islope=1, nslope 834 818 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)) 836 820 call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope)) 837 821 ENDDO 838 822 839 823 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) 841 825 842 826 !------------------------ … … 853 837 if (co2glaciersflow) then 854 838 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) 856 840 endif 857 841 … … 860 844 DO islope=1, nslope 861 845 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)) 863 847 call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) 864 848 ENDDO … … 881 865 DO ig = 1,ngrid 882 866 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 co2ice867 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 884 868 if(latitude_deg(ig).gt.0) then 885 869 do ig_loop=ig,ngrid 886 870 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) then871 if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then 888 872 tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop) 889 873 bool_sublim=1 … … 898 882 do ig_loop=ig,1,-1 899 883 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) then884 if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then 901 885 tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop) 902 886 bool_sublim=1 … … 909 893 enddo 910 894 endif 911 initial_co2_ice _slope(ig,islope)=0912 if (( co2ice_slope(ig,islope).lt.1e-10).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then913 albedo _slope(ig,1,islope) = albedo_h2o_frost914 albedo _slope(ig,2,islope) = albedo_h2o_frost895 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 915 899 else 916 albedo _slope(ig,1,islope) = albedodat(ig)917 albedo _slope(ig,2,islope) = albedodat(ig)918 emis s_slope(ig,islope) = emissiv900 albedo(ig,1,islope) = albedodat(ig) 901 albedo(ig,2,islope) = albedodat(ig) 902 emis(ig,islope) = emissiv 919 903 endif 920 elseif( ( co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2904 elseif( (qsurf(ig,igcm_co2,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2 921 905 ave=0. 922 906 do t=1,timelen 923 if(co2_ice_GCM _slope(ig,islope,t).gt.1e-3) then907 if(co2_ice_GCM(ig,islope,t).gt.1e-3) then 924 908 ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.)) 925 909 else … … 938 922 do ig = 1,ngrid 939 923 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)) 941 925 enddo 942 926 enddo … … 944 928 DO islope=1, nslope 945 929 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)) 947 931 ENDDO 948 932 … … 990 974 991 975 ! 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, & 993 977 porefillingice_depth,porefillingice_thickness,TI_PEM) 994 978 995 979 ! II_d.5 Update the mass of the regolith adsorbded 996 980 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,:), & 998 982 tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & 999 983 h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded) … … 1013 997 1014 998 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,& 1016 1000 global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope) 1017 1001 … … 1027 1011 1028 1012 !------------------------ 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) 1033 1017 1034 1018 year_iter=year_iter+dt_pem … … 1084 1068 1085 1069 ! III_a.1 Ice update (for startfi) 1070 #ifdef CPP_STD 1086 1071 ! Co2 ice 1087 1072 DO ig = 1,ngrid … … 1092 1077 cos(pi*def_slope_mean(islope)/180.) 1093 1078 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 1096 1081 #endif 1097 ENDDO ! of DO ig=1,ngrid1098 1082 1099 1083 ! H2O ice … … 1102 1086 watercap_sum=0. 1103 1087 DO islope=1,nslope 1104 watercap_s lope_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.))) then1106 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.)) 1107 1091 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. 1110 1094 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. 1113 1097 enddo 1114 1098 water_reservoir(ig)=water_reservoir(ig)+watercap_sum … … 1116 1100 enddo 1117 1101 1118 DO ig = 1,ngrid1119 qsurf(ig,igcm_h2o_ice) = 0.1120 DO islope = 1,nslope1121 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 ENDDO1125 ENDDO ! of DO ig=1,ngrid1126 1127 1102 DO ig=1,ngrid 1128 1103 DO islope=1,nslope 1129 if(qsurf _slope(ig,igcm_h2o_ice,islope).GT.500) then1104 if(qsurf(ig,igcm_h2o_ice,islope).GT.500) then 1130 1105 watercaptag(ig)=.true. 1131 qsurf _slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-2501106 qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)-250 1132 1107 water_reservoir(ig)=water_reservoir(ig)+250 1133 1108 endif … … 1136 1111 1137 1112 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 1155 1120 1156 1121 ! III_a.2 Tsoil update (for startfi) 1157 1122 1158 1123 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) 1162 1126 endif !soil_pem 1163 1164 #ifndef CPP_STD1165 DO ig = 1,ngrid1166 DO iloop = 1,nsoilmx1167 tsoil(ig,iloop) = 0.1168 inertiesoil(ig,iloop) = 0.1169 DO islope = 1,nslope1170 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 ENDDO1175 ENDDO1176 ENDDO ! of DO ig=1,ngrid1177 1178 ! III_a.3 Surface optical properties (for startfi)1179 1180 DO ig = 1,ngrid1181 DO l = 1,21182 albedo(ig,l) =0.1183 DO islope = 1,nslope1184 albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) &1185 *subslope_dist(ig,islope)1186 ENDDO1187 ENDDO1188 ENDDO1189 1190 DO ig = 1,ngrid1191 emis(ig) =0.1192 DO islope = 1,nslope1193 emis(ig)= emis(ig)+emiss_slope(ig,islope) &1194 *subslope_dist(ig,islope)1195 ENDDO1196 ENDDO1197 1198 DO ig = 1,ngrid1199 tsurf(ig) = 0.1200 DO islope = 1,nslope1201 tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)*tsurf_slope(ig,islope))**4 &1202 * subslope_dist(ig,islope)1203 ENDDO1204 tsurf(ig) = tsurf(ig)**(0.25)/emis(ig)1205 ENDDO ! of DO ig=1,ngrid1206 1207 #endif1208 1127 1209 1128 ! III_a.4 Pressure (for start) … … 1250 1169 DO ig=1,ngrid 1251 1170 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") ) then1171 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 1253 1172 extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1)) 1254 1173 q(ig,l,nnq)=1. … … 1295 1214 nsoilmx,ngrid,nlayer,nq, & 1296 1215 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, & 1299 1217 subslope_dist) 1300 1218 1301 1219 call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & 1302 1220 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) 1308 1224 #else 1309 1225 call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & … … 1346 1262 deallocate(q_co2_PEM_phys) 1347 1263 deallocate(q_h2o_PEM_phys) 1348 deallocate(co2_ice_GCM _slope)1264 deallocate(co2_ice_GCM) 1349 1265 deallocate(watersurf_density_ave) 1350 1266 deallocate(watersoil_density_timeseries) -
trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90
r2944 r2963 91 91 print *, "Downloading data for vmr co2..." 92 92 93 CALL get_var3("co2_ cropped" ,q_co2_dyn)93 CALL get_var3("co2_layer1" ,q_co2_dyn) 94 94 95 95 print *, "Downloading data for vmr co2 done" 96 96 print *, "Downloading data for vmr h20..." 97 97 98 CALL get_var3("h2o_ cropped" ,q_h2o_dyn)98 CALL get_var3("h2o_layer1" ,q_h2o_dyn) 99 99 100 100 print *, "Downloading data for vmr h2o done" … … 143 143 if(soil_pem) then 144 144 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 147 DO islope=1,nslope 148 write(num,fmt='(i2.2)') islope 149 call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:)) 150 ENDDO 151 152 print *, "Downloading data for soiltemp_slope done" 153 153 154 154 print *, "Downloading data for watersoil_density ..." … … 182 182 183 183 if(soil_pem) then 184 call get_var4(" tsoil",tsoil_gcm_dyn(:,:,:,1,:))184 call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:)) 185 185 endif !soil_pem 186 186 endif !nslope=1 -
trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
r2894 r2963 16 16 #ifndef CPP_STD 17 17 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 19 19 #else 20 20 use planete_mod, only: e_elips, obliquit, timeperi … … 45 45 real :: timeperi_ls ! time of peri in ls 46 46 integer ilask ! Loop variable 47 real :: halfaxe ! Million of km 48 real :: unitastr 47 49 48 50 … … 73 75 endif 74 76 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 75 84 timeperi_ls=lsplask(ilask+1)+(lsplask(ilask)-lsplask(ilask+1))*(Year-yearlask(ilask+1))/(yearlask(ilask)-yearlask(ilask+1)) 85 timeperi_ls=modulo(timeperi_ls,360.) 76 86 endif 77 87 exit 78 88 endif 79 89 enddo 90 halfaxe=227.94 80 91 timeperi=timeperi_ls*2*pi/360 92 periheli = halfaxe*(1-e_elips) 93 aphelie = halfaxe*(1+e_elips) 94 unitastr=149.597927 95 p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr 96 call call_dayperi(timeperi,e_elips,peri_day,year_day) 81 97 82 98 print *, "recomp_orb_param, Final year of the PEM run:", Year -
trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90
r2855 r2963 51 51 integer :: numyear 52 52 53 54 55 56 57 53 DO numyear=1, 2 58 54 write(*,*) 'numyear',numyear … … 61 57 !integer :: ncid, status 62 58 !... 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)59 status = nf90_open(path = "data2reshape"//str2//".nc", mode = nf90_nowrite, ncid = ncid1) 60 if(status /= nf90_noerr) call handle_err(status) 61 62 status = nf90_create(path = "datareshaped"//str2//".nc", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2) 67 63 if(status /= nf90_noerr) call handle_err(status) 68 64 -
trunk/LMDZ.MARS/deftank/context_lmdz_physics.xml
r2934 r2963 129 129 130 130 <!-- diurnal averages outputs; output_freq is every day --> 131 <file id="diurnalave "132 name="Xdiurnalave "131 <file id="diurnalave_" 132 name="Xdiurnalave_" 133 133 output_freq="1d" 134 134 time_units="days" … … 157 157 </field_group> 158 158 </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 159 245 </file_definition> 160 246 -
trunk/LMDZ.MARS/deftank/field_def_physics_mars.xml
r2934 r2963 68 68 long_name="Surface emissivity of slope 01" 69 69 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="" /> 70 88 <field id="albedo" 71 89 long_name="Albedo of the surface" … … 74 92 long_name="Albedo of the surface for slope 01" 75 93 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="" /> 77 112 <field id="local_time" 78 113 long_name="Local time" … … 88 123 long_name="Surface Temperature of slope 01" 89 124 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" /> 90 143 <field id="temp_layer1" 91 144 long_name="Temperature in layer 1" … … 102 155 long_name="Longwave radiation at the surface on slope 01" 103 156 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" /> 104 175 <field id="fluxtop_lw" 105 176 long_name="Longwave radiation at the top of the atmosphere" … … 111 182 <field id="fluxsurf_dn_sw_slope01" 112 183 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" 113 202 unit="W.m-2" /> 114 203 <field id="fluxtop_up_sw" … … 145 234 long_name="tendency due to vertical diffusion of background dust on surface" 146 235 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" /> 147 254 <field id="dqsdifrdsq" 148 255 long_name="tendency due to vertical diffusion of stormdust on surface" 149 256 unit="kg.m-2.s-1" /> 150 257 <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" 151 276 long_name="tendency due to vertical diffusion of stormdust on surface" 152 277 unit="kg.m-2.s-1" /> … … 185 310 <field id="watercap_slope01" 186 311 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 188 332 <field id="surf_h2o_lh" 189 333 long_name="Ground ice latent heat flux" … … 213 357 long_name="Mass of water ice on the surface of slope 01" 214 358 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 215 406 216 407 <!-- CO2 cycle --> … … 221 412 long_name="CO2 ice thickness on slope 01" 222 413 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 223 433 224 434 <field id="co2condens_zfallice" … … 463 673 long_name="Soil temperature for slope 01" 464 674 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" /> 465 693 <field id="inertiedat" 466 694 long_name="Soil thermal inertia" 467 695 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 468 739 </field_group> 469 740 -
trunk/LMDZ.MARS/libf/phymars/call_dayperi.F
r2448 r2963 17 17 c Lsperi solar longitude (Ls) of perohelion (rad) 18 18 c e_elips Excentricity 19 c real year_day ! number of sols per Mars yar 19 20 c 20 21 c output 21 22 c ------ 22 23 c dayperi Martian date at perihelion (sol) 23 c real year_day ! number of sols per Mars yar 24 24 25 c----------------------------------------------------------------------- 25 26 -
trunk/LMDZ.MARS/libf/phymars/iniorbit.F
r1382 r2963 68 68 69 69 timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) 70 if(timeperi.lt.0) timeperi=-timeperi 70 71 PRINT*,'iniorbit: Perihelion solar long. Ls (deg)=', 71 72 & 360.-timeperi*180./pi -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r2953 r2963 545 545 logical :: write_restart 546 546 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 547 558 c======================================================================= 548 559 pdq(:,:,:) = 0. … … 789 800 call update_xios_timestep 790 801 #endif 802 803 791 804 792 805 … … 1498 1511 dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:) 1499 1512 ENDDO 1513 1500 1514 call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif, 1501 1515 & albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp) … … 1559 1573 ENDDO 1560 1574 ENDDO 1575 1561 1576 IF (turb_resolved) THEN 1562 1577 write(*,*) 'Turbulent-resolving mode !' … … 2160 2175 endif ! of if (callthermos) 2161 2176 #endif 2177 2162 2178 c----------------------------------------------------------------------- 2163 2179 c 11. Carbon dioxide condensation-sublimation: … … 3822 3838 endif !not.scavenging 3823 3839 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 3824 3880 !PREVIOUSLY IN 1D ONLY 3825 3881 -
trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90
r2953 r2963 532 532 'lapse rate in the rocket dust storm', & 533 533 '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(:,:)) 537 537 ! call write_output('scheme','which scheme',& 538 538 ! ' ',real(scheme(:))) -
trunk/LMDZ.MARS/util/startarchive2icosa/rearrange_startphy.f90
r2941 r2963 270 270 endif 271 271 enddo 272 273 print *, "All is well that ends well" 272 274 273 275 END PROGRAM rearrange_startphy
Note: See TracChangeset
for help on using the changeset viewer.