Changeset 3143
- Timestamp:
- Nov 29, 2023, 9:42:08 AM (12 months ago)
- Location:
- trunk
- Files:
-
- 8 edited
-
LMDZ.COMMON/libf/evolution/changelog.txt (modified) (1 diff)
-
LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90 (modified) (1 diff)
-
LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90 (modified) (1 diff)
-
LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90 (modified) (4 diffs)
-
LMDZ.COMMON/libf/evolution/pem.F90 (modified) (32 diffs)
-
LMDZ.COMMON/libf/evolution/pemetat0.F90 (modified) (3 diffs)
-
LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 (modified) (4 diffs)
-
LMDZ.MARS/deftank/pem/concat_diagpem.sh (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3136 r3143 145 145 == 22/11/2023 == JBC 146 146 Update of files in the deftank: addition of variables in the xml definition files, inclusion of "callphys.def" into "run_PEM.def" and minor typo in "launch_pem.sh". 147 148 == 29/11/2023 == JBC 149 'Watercap' has been removed from the water ice evolution since 'water_reservoir' does already the job + Some cleanings to simplify the code. 150 151 /!\ Commit for the PEM management of h2o ice before a rework of ice management in the PEM! -
trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90
r3123 r3143 1 modulecomsoil_h_PEM1 MODULE comsoil_h_PEM 2 2 3 3 implicit none 4 integer, parameter :: nsoilmx_PEM = 68 ! number of layers in the PEM5 real,save,allocatable,dimension(:) :: layer_PEM ! soil layer depths [m]6 real,save,allocatable,dimension(:) :: mlayer_PEM ! soil mid-layer depths [m]7 real,save,allocatable,dimension(:,:,:) :: TI_PEM ! soil thermal inertia [SI]8 real,save,allocatable,dimension(:,:) :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI]9 ! variables (FC: built in firstcall in soil.F)10 REAL,SAVE,ALLOCATABLE :: tsoil_PEM(:,:,:) ! sub-surface temperatures [K]11 real,save,allocatable :: mthermdiff_PEM(:,:) ! (FC) mid-layer thermal diffusivity [SI]12 real,save,allocatable :: thermdiff_PEM(:,:) ! (FC) inter-layer thermal diffusivity [SI]13 real,save,allocatable :: coefq_PEM(:) ! (FC) q_{k+1/2} coefficients [SI]14 real,save,allocatable :: coefd_PEM(:,:) ! (FC) d_k coefficients [SI]15 real,save,allocatable :: alph_PEM(:,:) ! (FC) alpha_k coefficients [SI]16 real,save,allocatable :: beta_PEM(:,:) ! beta_k coefficients [SI]17 real,save :: mu_PEM ! mu coefficient [SI]18 real,save :: fluxgeo ! Geothermal flux [W/m^2]19 real,save :: depth_breccia ! Depth at which we have breccia [m]20 real,save :: depth_bedrock ! Depth at which we have bedrock [m]21 integer,save :: index_breccia ! last index of the depth grid before having breccia22 integer,save :: index_bedrock ! last index of the depth grid before having bedrock23 LOGICAL soil_pem ! True by default, to run with the subsurface physic. Read in pem.def24 real,save,allocatable,dimension(:) :: water_reservoir ! Large reserve of water [kg/m^2]25 real,save :: water_reservoir_nom ! Nominal value for the large reserve of water [kg/m^2]26 logical, save :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure27 4 5 integer, parameter :: nsoilmx_PEM = 68 ! number of layers in the PEM 6 real, allocatable, dimension(:), save :: layer_PEM ! soil layer depths [m] 7 real, allocatable, dimension(:), save :: mlayer_PEM ! soil mid-layer depths [m] 8 real, allocatable, dimension(:,:,:), save :: TI_PEM ! soil thermal inertia [SI] 9 real, allocatable, dimension(:,:), save :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI] 10 ! Variables (FC: built in firstcall in soil.F) 11 real, allocatable, dimension(:,:,:), save :: tsoil_PEM ! sub-surface temperatures [K] 12 real, allocatable, dimension(:,:), save :: mthermdiff_PEM ! (FC) mid-layer thermal diffusivity [SI] 13 real, allocatable, dimension(:,:), save :: thermdiff_PEM ! (FC) inter-layer thermal diffusivity [SI] 14 real, allocatable, dimension(:), save :: coefq_PEM ! (FC) q_{k+1/2} coefficients [SI] 15 real, allocatable, dimension(:,:), save :: coefd_PEM ! (FC) d_k coefficients [SI] 16 real, allocatable, dimension(:,:), save :: alph_PEM ! (FC) alpha_k coefficients [SI] 17 real, allocatable, dimension(:,:), save :: beta_PEM ! beta_k coefficients [SI] 18 real, save :: mu_PEM ! mu coefficient [SI] 19 real, save :: fluxgeo ! Geothermal flux [W/m^2] 20 real, save :: depth_breccia ! Depth at which we have breccia [m] 21 real, save :: depth_bedrock ! Depth at which we have bedrock [m] 22 integer, save :: index_breccia ! last index of the depth grid before having breccia 23 integer, save :: index_bedrock ! last index of the depth grid before having bedrock 24 logical :: soil_pem ! True by default, to run with the subsurface physic. Read in pem.def 25 real, allocatable, dimension(:), save :: water_reservoir ! Large reserve of water [kg/m^2] 26 real, save :: water_reservoir_nom ! Nominal value for the large reserve of water [kg/m^2] 27 logical, save :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure 28 29 !======================================================================= 28 30 contains 31 !======================================================================= 29 32 30 subroutine ini_comsoil_h_PEM(ngrid,nslope) 31 32 implicit none 33 integer,intent(in) :: ngrid ! number of atmospheric columns 34 integer,intent(in) :: nslope ! number of slope within a mesh 33 SUBROUTINE ini_comsoil_h_PEM(ngrid,nslope) 35 34 36 allocate(layer_PEM(nsoilmx_PEM)) 37 allocate(mlayer_PEM(0:nsoilmx_PEM-1)) 38 allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope)) 39 allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope)) 40 allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1)) 41 allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1)) 42 allocate(coefq_PEM(0:nsoilmx_PEM-1)) 43 allocate(coefd_PEM(ngrid,nsoilmx_PEM-1)) 44 allocate(alph_PEM(ngrid,nsoilmx_PEM-1)) 45 allocate(beta_PEM(ngrid,nsoilmx_PEM-1)) 46 allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) 47 allocate(water_reservoir(ngrid)) 48 end subroutine ini_comsoil_h_PEM 35 implicit none 49 36 37 integer, intent(in) :: ngrid ! number of atmospheric columns 38 integer, intent(in) :: nslope ! number of slope within a mesh 50 39 51 subroutine end_comsoil_h_PEM 40 allocate(layer_PEM(nsoilmx_PEM)) 41 allocate(mlayer_PEM(0:nsoilmx_PEM-1)) 42 allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope)) 43 allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope)) 44 allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1)) 45 allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1)) 46 allocate(coefq_PEM(0:nsoilmx_PEM-1)) 47 allocate(coefd_PEM(ngrid,nsoilmx_PEM-1)) 48 allocate(alph_PEM(ngrid,nsoilmx_PEM-1)) 49 allocate(beta_PEM(ngrid,nsoilmx_PEM-1)) 50 allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) 51 allocate(water_reservoir(ngrid)) 52 52 53 implicit none 53 END SUBROUTINE ini_comsoil_h_PEM 54 54 55 if (allocated(layer_PEM)) deallocate(layer_PEM) 56 if (allocated(mlayer_PEM)) deallocate(mlayer_PEM) 57 if (allocated(TI_PEM)) deallocate(TI_PEM) 58 if (allocated(tsoil_PEM)) deallocate(tsoil_PEM) 59 if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM) 60 if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM) 61 if (allocated(coefq_PEM)) deallocate(coefq_PEM) 62 if (allocated(coefd_PEM)) deallocate(coefd_PEM) 63 if (allocated(alph_PEM)) deallocate(alph_PEM) 64 if (allocated(beta_PEM)) deallocate(beta_PEM) 65 if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM) 66 if (allocated(water_reservoir)) deallocate(water_reservoir) 67 end subroutine end_comsoil_h_PEM 55 !======================================================================= 68 56 69 end module comsoil_h_PEM 57 SUBROUTINE end_comsoil_h_PEM 58 59 implicit none 60 61 if (allocated(layer_PEM)) deallocate(layer_PEM) 62 if (allocated(mlayer_PEM)) deallocate(mlayer_PEM) 63 if (allocated(TI_PEM)) deallocate(TI_PEM) 64 if (allocated(tsoil_PEM)) deallocate(tsoil_PEM) 65 if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM) 66 if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM) 67 if (allocated(coefq_PEM)) deallocate(coefq_PEM) 68 if (allocated(coefd_PEM)) deallocate(coefd_PEM) 69 if (allocated(alph_PEM)) deallocate(alph_PEM) 70 if (allocated(beta_PEM)) deallocate(beta_PEM) 71 if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM) 72 if (allocated(water_reservoir)) deallocate(water_reservoir) 73 74 END SUBROUTINE end_comsoil_h_PEM 75 76 END MODULE comsoil_h_PEM -
trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90
r3130 r3143 38 38 39 39 ! Conversion H2O/CO2 frost to perennial frost and vice versa 40 real, parameter :: threshold_ water_frost2perennial = 1000.!~ 1 m41 real, parameter :: threshold_co2_frost2perennial = 16000. 40 real, parameter :: threshold_h2o_frost2perennial = 1000. !~ 1 m 41 real, parameter :: threshold_co2_frost2perennial = 16000. !~ 10 m 42 42 43 43 END MODULE constants_marspem_mod 44 -
trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90
r3130 r3143 3 3 implicit none 4 4 5 ! ---------------------------------------------------------------------- 5 6 contains 7 ! ---------------------------------------------------------------------- 6 8 7 9 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 44 46 STOPPING = .false. 45 47 46 ! Computation of the present surface of water ice sublimating48 ! Computation of the present surface of h2o ice 47 49 present_surf = 0. 48 50 do i = 1,ngrid 49 51 do islope = 1,nslope 50 !if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope) 51 if (initial_h2o_ice(i,islope) > 0.5) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope) 52 if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope) 52 53 enddo 53 54 enddo 54 55 55 56 ! Check of the criterion 56 if (present_surf < ini_surf*(1. - water_ice_criterion) .or. present_surf > ini_surf*(1. + water_ice_criterion)) then57 if (present_surf < ini_surf*(1. - water_ice_criterion)) then 57 58 STOPPING = .true. 58 59 write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold" 60 write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", present_surf < ini_surf*(1. - water_ice_criterion) 59 61 write(*,*) "Current surface of water ice sublimating =", present_surf 60 62 write(*,*) "Initial surface of water ice sublimating =", ini_surf 61 63 write(*,*) "Percentage of change accepted =", water_ice_criterion*100 62 write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", (present_surf < ini_surf*(1. - water_ice_criterion)) 64 endif 65 if (present_surf > ini_surf*(1. + water_ice_criterion)) then 66 STOPPING = .true. 67 write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold" 68 write(*,*) "present_surf > ini_surf*(1. + water_ice_criterion)", present_surf > ini_surf*(1. + water_ice_criterion) 69 write(*,*) "Current surface of water ice sublimating =", present_surf 70 write(*,*) "Initial surface of water ice sublimating =", ini_surf 71 write(*,*) "Percentage of change accepted =", water_ice_criterion*100 63 72 endif 64 73 … … 105 114 STOPPING_ps = .false. 106 115 107 ! Computation of the actual surface116 ! Computation of the present surface of co2 ice 108 117 present_surf = 0. 109 118 do i = 1,ngrid … … 114 123 115 124 ! Check of the criterion 116 if (present_surf < ini_surf*(1. - co2_ice_criterion) .or. present_surf > ini_surf*(1. + co2_ice_criterion)) then125 if (present_surf < ini_surf*(1. - co2_ice_criterion)) then 117 126 STOPPING_ice = .true. 118 write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold" 127 write(*,*) "Reason of stopping: the surface of co2 ice sublimating reached the threshold" 128 write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", present_surf < ini_surf*(1. - co2_ice_criterion) 119 129 write(*,*) "Current surface of co2 ice sublimating =", present_surf 120 130 write(*,*) "Initial surface of co2 ice sublimating =", ini_surf 121 131 write(*,*) "Percentage of change accepted =", co2_ice_criterion*100. 122 write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", (present_surf < ini_surf*(1. - co2_ice_criterion)) 132 endif 133 if (present_surf > ini_surf*(1. + co2_ice_criterion)) then 134 STOPPING_ice = .true. 135 write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold" 136 write(*,*) "present_surf > ini_surf*(1. + co2_ice_criterion)", present_surf > ini_surf*(1. + co2_ice_criterion) 137 write(*,*) "Current surface of co2 ice sublimating =", present_surf 138 write(*,*) "Initial surface of co2 ice sublimating =", ini_surf 139 write(*,*) "Percentage of change accepted =", co2_ice_criterion*100. 123 140 endif 124 141 125 142 if (abs(ini_surf) < 1.e-5) STOPPING_ice = .false. 126 143 127 if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion) .or. global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then144 if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)) then 128 145 STOPPING_ps = .true. 129 146 write(*,*) "Reason of stopping: the global pressure reach the threshold" 147 write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion) 130 148 write(*,*) "Current global pressure =", global_ave_press_new 131 149 write(*,*) "PCM global pressure =", global_ave_press_GCM 132 150 write(*,*) "Percentage of change accepted =", ps_criterion*100. 133 write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)) 151 endif 152 if (global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then 153 STOPPING_ps = .true. 154 write(*,*) "Reason of stopping: the global pressure reach the threshold" 155 write(*,*) "global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)", global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion) 156 write(*,*) "Current global pressure =", global_ave_press_new 157 write(*,*) "PCM global pressure =", global_ave_press_GCM 158 write(*,*) "Percentage of change accepted =", ps_criterion*100. 134 159 endif 135 160 -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3130 r3143 1 1 !------------------------ 2 2 ! I Initialization 3 ! I_a R EADrun.def4 ! I_b R EADof start_evol.nc and starfi_evol.nc3 ! I_a Read run.def 4 ! I_b Read of start_evol.nc and starfi_evol.nc 5 5 ! I_c Subslope parametrisation 6 ! I_d R EADPCM data and convert to the physical grid6 ! I_d Read PCM data and convert to the physical grid 7 7 ! I_e Initialization of the PEM variable and soil 8 ! I_f Compute tendencies & Save initial situation8 ! I_f Compute tendencies 9 9 ! I_g Save initial PCM situation 10 10 ! I_h Read the startpem.nc … … 42 42 use criterion_pem_stop_mod, only: criterion_waterice_stop, criterion_co2_stop 43 43 use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, & 44 m_noco2, threshold_ water_frost2perennial, threshold_co2_frost2perennial44 m_noco2, threshold_h2o_frost2perennial, threshold_co2_frost2perennial 45 45 use evol_co2_ice_s_mod, only: evol_co2_ice_s 46 46 use evol_h2o_ice_s_mod, only: evol_h2o_ice_s … … 86 86 use paleoclimate_mod, only: albedo_perennialco2 87 87 use comcstfi_h, only: pi, rad, g, cpp, mugaz, r 88 use surfini_mod, only: surfini 88 89 #else 89 90 use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names … … 131 132 132 133 ! Variables to read start.nc 133 character( len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsingthe PEM134 character(*), parameter :: start_name = "start_evol.nc" ! Name of the file used to initialize the PEM 134 135 135 136 ! Dynamic variables … … 138 139 real, dimension(ip1jmp1,llm) :: teta ! temperature potentielle 139 140 real, dimension(:,:,:), allocatable :: q ! champs advectes 140 real, dimension(ip1jmp1) :: ps ! pression 141 real, dimension(ip1jmp1) :: ps ! pression au sol 141 142 real, dimension(:), allocatable :: ps_start_GCM ! (ngrid) pression au sol 142 real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) ! pression 143 real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées 143 144 real, dimension(ip1jmp1,llm) :: masse ! masse d'air 144 145 real, dimension(ip1jmp1) :: phis ! geopotentiel au sol … … 146 147 147 148 ! Variables to read starfi.nc 148 character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" ! Name of the file used for initialsingthe PEM149 character(2) 150 integer 151 integer 152 integer 153 integer 149 character(*), parameter :: startfi_name = "startfi_evol.nc" ! Name of the file used to initialize the PEM 150 character(2) :: str2 151 integer :: ncid, varid, status ! Variable for handling opening of files 152 integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files 153 integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files 154 integer :: apvarid, bpvarid ! Variable ID for Netcdf files 154 155 155 156 ! Variables to read starfi.nc and write restartfi.nc 156 real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAMEand written in restartfi157 real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAMEand written in restartfi158 real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAMEand written in restartfi157 real, dimension(:), allocatable :: longitude ! Longitude read in startfi_name and written in restartfi 158 real, dimension(:), allocatable :: latitude ! Latitude read in startfi_name and written in restartfi 159 real, dimension(:), allocatable :: cell_area ! Cell_area read in startfi_name and written in restartfi 159 160 real :: Total_surface ! Total surface of the planet 160 161 … … 187 188 real, dimension(:,:), allocatable :: co2frost_ini ! Initial amount of co2 frost (at the start of the PEM run) 188 189 real, dimension(:,:), allocatable :: perennial_co2ice_ini ! Initial amoun of perennial co2 ice (at the start of the PEM run) 189 190 191 190 192 191 ! Variables for slopes … … 243 242 real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice 244 243 real :: albedo_h2o_frost ! albedo of h2o frost 245 real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimensi ion when reading form generic246 real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimensi ion when reading form generic247 real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimensi ion when reading form generic248 real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimensi ion when reading form generic249 real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimensi ion when reading form generic244 real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic 245 real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic 246 real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic 247 real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic 248 real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic 250 249 real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case 251 250 real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case … … 259 258 260 259 #ifdef CPP_1D 261 integer 262 integer, parameter 263 integer 260 integer :: nsplit_phys 261 integer, parameter :: jjm_value = jjm - 1 262 integer :: ierr 264 263 265 264 ! Dummy variables to use the subroutine 'init_testphys1d' … … 272 271 real, dimension(nlayer + 1) :: plev 273 272 #else 274 integer, parameter 275 real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_startand written in restart276 real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_startand written in restart273 integer, parameter :: jjm_value = jjm 274 real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart 275 real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart 277 276 #endif 278 277 … … 302 301 !------------------------ 303 302 ! I Initialization 304 ! I_a R EADrun.def303 ! I_a Read run.def 305 304 !------------------------ 306 305 #ifndef CPP_1D … … 315 314 allocate(longitude(1),latitude(1),cell_area(1)) 316 315 317 therestart1D = .false. 316 therestart1D = .false. ! Default value 318 317 inquire(file = 'start1D_evol.txt',exist = therestart1D) 319 318 if (.not. therestart1D) then … … 321 320 error stop 'Initialization cannot be done for the 1D PEM.' 322 321 endif 323 therestartfi = .false. 322 therestartfi = .false. ! Default value 324 323 inquire(file = 'startfi_evol.nc',exist = therestartfi) 325 324 if (.not. therestartfi) then … … 339 338 !------------------------ 340 339 ! I Initialization 341 ! I_b R EADof start_evol.nc and starfi_evol.nc342 !------------------------ 343 ! I_b.1 R EADstart_evol.nc340 ! I_b Read of start_evol.nc and starfi_evol.nc 341 !------------------------ 342 ! I_b.1 Read start_evol.nc 344 343 allocate(ps_start_GCM(ngrid)) 345 344 #ifndef CPP_1D 346 call dynetat0( FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0)345 call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0) 347 346 348 347 call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM) … … 353 352 allocate(ap(nlayer + 1)) 354 353 allocate(bp(nlayer + 1)) 355 status = nf90_open( FILE_NAME_start,NF90_NOWRITE,ncid)354 status = nf90_open(start_name,NF90_NOWRITE,ncid) 356 355 status = nf90_inq_varid(ncid,"ap",apvarid) 357 356 status = nf90_get_var(ncid,apvarid,ap) … … 368 367 ! In the PCM, these values are given to the physic by the dynamic. 369 368 ! Here we simply read them in the startfi_evol.nc file 370 status = nf90_open( FILE_NAME, NF90_NOWRITE, ncid)369 status = nf90_open(startfi_name, NF90_NOWRITE, ncid) 371 370 372 371 status = nf90_inq_varid(ncid,"longitude",lonvarid) … … 384 383 status = nf90_close(ncid) 385 384 386 ! I_b.2 R EADstartfi_evol.nc385 ! I_b.2 Read startfi_evol.nc 387 386 ! First we read the initial state (starfi.nc) 388 387 #ifndef CPP_STD 389 call phyetat0( FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &388 call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, & 390 389 tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, & 391 390 watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist) … … 394 393 where (qsurf < 0.) qsurf = 0. 395 394 396 call surfini(ngrid, qsurf)395 call surfini(ngrid,nslope,qsurf) 397 396 #else 398 397 call phys_state_var_init(nq) … … 414 413 allocate(albedo(ngrid,2,1)) 415 414 allocate(inertiesoil(ngrid,nsoilmx,1)) 416 call phyetat0(.true.,ngrid,nlayer, FILE_NAME,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &415 call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, & 417 416 tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, & 418 417 qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, & … … 449 448 if (noms(nnq) == "h2o_vap") then 450 449 igcm_h2o_vap = nnq 451 mmol(igcm_h2o_vap) =18.450 mmol(igcm_h2o_vap) = 18. 452 451 endif 453 452 if (noms(nnq) == "co2") igcm_co2 = nnq … … 480 479 !------------------------ 481 480 ! I Initialization 482 ! I_d R EADPCM data and convert to the physical grid481 ! I_d Read PCM data and convert to the physical grid 483 482 !------------------------ 484 483 ! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value … … 552 551 !------------------------ 553 552 ! I Initialization 554 ! I_f Compute tendencies & Save initial situation553 ! I_f Compute tendencies 555 554 !------------------------ 556 555 allocate(tendencies_co2_ice(ngrid,nslope)) … … 603 602 write(*,*) "Total initial surface of co2 ice sublimating (slope) =", ini_surf_co2 604 603 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o 605 write(*,*) "Total surface of the planet =", Total_surface604 write(*,*) "Total surface of the planet =", Total_surface 606 605 allocate(zplev_gcm(ngrid,nlayer + 1)) 607 606 … … 617 616 global_ave_press_GCM = global_ave_press_old 618 617 global_ave_press_new = global_ave_press_old 619 write(*,*) "Initial global average pressure =", global_ave_press_GCM618 write(*,*) "Initial global average pressure =", global_ave_press_GCM 620 619 621 620 !------------------------ … … 638 637 do ig = 1,ngrid 639 638 do islope = 1,nslope 640 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water cap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)639 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) 641 640 qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) + perennial_co2ice(ig,islope) 642 641 enddo … … 657 656 enddo 658 657 659 write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded660 write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded658 write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded 659 write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded 661 660 endif ! adsorption 662 661 deallocate(tsurf_ave_yr1) … … 735 734 q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ & 736 735 (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) 737 if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30 738 if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1. 736 if (q_h2o_PEM_phys(ig,t) < 0) then 737 q_h2o_PEM_phys(ig,t) = 1.e-30 738 else if (q_h2o_PEM_phys(ig,t) > 1) then 739 q_h2o_PEM_phys(ig,t) = 1. 740 endif 739 741 enddo 740 742 enddo … … 749 751 if (q_co2_PEM_phys(ig,t) < 0) then 750 752 q_co2_PEM_phys(ig,t) = 1.e-30 751 else if (q_co2_PEM_phys(ig,t) > 1) then753 else if (q_co2_PEM_phys(ig,t) > 1) then 752 754 q_co2_PEM_phys(ig,t) = 1. 753 755 endif … … 898 900 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) 899 901 900 ! II_d.5 Update the mass of the regolith adsorb ded902 ! II_d.5 Update the mass of the regolith adsorbed 901 903 if (adsorption_pem) then 902 904 call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, & … … 985 987 exit 986 988 else 987 write(*,*) " Wecontinue!"989 write(*,*) "The PEM can continue!" 988 990 write(*,*) "Number of iterations of the PEM: year_iter =", year_iter 989 991 write(*,*) "Number of simulated Martian years: i_myear =", i_myear … … 1007 1009 watercap_sum = 0. 1008 1010 do islope = 1,nslope 1009 if (qsurf(ig,igcm_h2o_ice,islope) > (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy 1010 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.)) ! put back ancien frost 1011 ! water_reservoir and watercap have not changed since PCM call: here we check if we have accumulate frost or not. 1012 if (qsurf(ig,igcm_h2o_ice,islope) > water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) then 1013 ! 1st case: more ice than initialy. Put back ancien frost 1014 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) 1011 1015 else 1012 ! 2nd case: we have sublimate ice: then let's put qsurf = 0.and add the difference in watercap1013 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.))1016 ! 2nd case: ice sublimated. Then let's put qsurf = 0 and add the difference in watercap 1017 watercap(ig,islope) = qsurf(ig,igcm_h2o_ice,islope) 1014 1018 qsurf(ig,igcm_h2o_ice,islope) = 0. 1015 1019 endif … … 1027 1031 enddo 1028 1032 if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. 1029 if (water_sum > threshold_ water_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.11033 if (water_sum > threshold_h2o_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1 1030 1034 watercaptag(ig) = .true. 1031 water_reservoir(ig) = water_reservoir(ig) + threshold_ water_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost1035 water_reservoir(ig) = water_reservoir(ig) + threshold_h2o_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost 1032 1036 do islope = 1,nslope 1033 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_ water_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.)1037 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_h2o_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.) 1034 1038 enddo 1035 1039 endif 1036 1040 else ! let's check that the infinite source of water disapear 1037 if ((water_sum + water_reservoir(ig)) < threshold_ water_frost2perennial) then1041 if ((water_sum + water_reservoir(ig)) < threshold_h2o_frost2perennial) then 1038 1042 watercaptag(ig) = .false. 1039 1043 do islope = 1,nslope -
trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
r3123 r3143 104 104 !0. Check if the start_PEM exist. 105 105 106 inquire(file =filename,exist =startpem_file)106 inquire(file = filename,exist = startpem_file) 107 107 108 108 write(*,*)'Is start PEM?',startpem_file … … 312 312 !. 5 water reservoir 313 313 #ifndef CPP_STD 314 water_reservoir = 0. 314 315 call get_field("water_reservoir",water_reservoir,found) 315 if(.not.found) then 316 write(*,*)'Pemetat0: failed loading <water_reservoir>' 317 write(*,*)'will reconstruct the values from watercaptag' 318 do ig=1,ngrid 319 if(watercaptag(ig)) then 320 water_reservoir(ig)=water_reservoir_nom 321 else 322 water_reservoir(ig)=0. 323 endif 324 enddo 316 if (.not. found) then 317 write(*,*)'Pemetat0: failed loading <water_reservoir>' 318 write(*,*)'will reconstruct the values from watercaptag' 319 where (watercaptag) water_reservoir = water_reservoir_nom 325 320 endif 326 321 #endif … … 481 476 !. e) water reservoir 482 477 #ifndef CPP_STD 483 do ig = 1,ngrid 484 if (watercaptag(ig)) then 485 water_reservoir(ig) = water_reservoir_nom 486 else 487 water_reservoir(ig) = 0. 488 endif 489 enddo 478 water_reservoir = 0. 479 where (watercaptag) water_reservoir = water_reservoir_nom 490 480 #endif 491 481 -
trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90
r3130 r3143 104 104 write(*,*) "Data for h2o_ice_s downloaded!" 105 105 106 call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:)) 107 write(*,*) "Data for tsurf downloaded!" 108 106 109 #ifndef CPP_STD 107 110 call get_var3("watercap",watercap(:,:,1,:)) … … 110 113 call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:)) 111 114 write(*,*) "Data for perennial_co2ice downloaded!" 112 #endif 113 114 call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:)) 115 write(*,*) "Data for tsurf downloaded!" 116 117 #ifndef CPP_STD 115 118 116 if (soil_pem) then 119 117 call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:)) … … 140 138 write(*,*) "Data for h2o_ice_s downloaded!" 141 139 140 do islope = 1,nslope 141 write(num,'(i2.2)') islope 142 call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:)) 143 enddo 144 write(*,*) "Data for tsurf downloaded!" 145 142 146 #ifndef CPP_STD 143 147 do islope = 1,nslope … … 152 156 enddo 153 157 write(*,*) "Data for perennial_co2ice downloaded!" 154 #endif 155 156 do islope = 1,nslope 157 write(num,'(i2.2)') islope 158 call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:)) 159 enddo 160 write(*,*) "Data for tsurf downloaded!" 161 162 #ifndef CPP_STD 158 163 159 if (soil_pem) then 164 160 do islope = 1,nslope -
trunk/LMDZ.MARS/deftank/pem/concat_diagpem.sh
r3088 r3143 58 58 done 59 59 60 echo "Concatenation of \"$directory/diagpem*.nc\" files into \"$ directory/$output_file\" is complete!"60 echo "Concatenation of \"$directory/diagpem*.nc\" files into \"$output_file\" is complete!"
Note: See TracChangeset
for help on using the changeset viewer.
