- Timestamp:
- Oct 2, 2023, 2:30:47 PM (21 months ago)
- Location:
- trunk
- Files:
-
- 1 deleted
- 6 edited
-
LMDZ.COMMON/libf/evolution/changelog.txt (modified) (1 diff)
-
LMDZ.COMMON/libf/evolution/init_pem1d_mod.F90 (deleted)
-
LMDZ.COMMON/libf/evolution/pem.F90 (modified) (31 diffs)
-
LMDZ.MARS/deftank/pem/launch_pem.sh (modified) (5 diffs)
-
LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 (modified) (11 diffs)
-
LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 (modified) (4 diffs)
-
LMDZ.MARS/libf/phymars/dyn1d/writerestart1D.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3063 r3065 77 77 == 26/09/2023 == JBC 78 78 Minor changes concerning the form of the code in the PEM. 79 Addition of a file "changelog.txt" in LMDZ.COMMON/libf/evolution/ specific to the PEM rather than using the one for Mars. Completed with changesets since 01/07/2023. 79 80 80 == 28/09/2023 == JBC81 Addition of a file "changelog.txt" in LMDZ.COMMON/libf/evolution/ specific to the PEM rather than using the one for Mars. Completed with changesets since 01/07/2023.81 == 02/10/2023 == JBC 82 Initialization of the PEM in 1D through the subroutine "init_testphys1d_mod.F90" + Some adaptations of the Mars PCM in 1D + Update of "launch_pem.sh" in deftank. -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3050 r3065 96 96 use physics_distribution_mod, only: init_physics_distribution 97 97 use mod_grid_phy_lmdz, only: regular_lonlat 98 use init_pem1d_mod, only: init_pem1d 98 use init_testphys1d_mod, only: init_testphys1d 99 use comvert_mod, only: ap, bp 99 100 #endif 100 101 … … 111 112 112 113 ! Same variable names as in the GCM 113 integer :: ngrid ! Number of physical grid points114 integer :: nlayer ! Number of vertical layer115 integer :: nq! Number of tracer116 integer :: day_ini! First day of the simulation117 real :: pday! Physical day118 real :: time_phys! Same as GCM119 real :: ptimestep! Same as GCM120 real :: ztime_fin! Same as GCM114 integer, parameter :: nlayer = llm ! Number of vertical layer 115 integer :: ngrid ! Number of physical grid points 116 integer :: nq ! Number of tracer 117 integer :: day_ini ! First day of the simulation 118 real :: pday ! Physical day 119 real :: time_phys ! Same as GCM 120 real :: ptimestep ! Same as GCM 121 real :: ztime_fin ! Same as GCM 121 122 122 123 ! Variables to read start.nc … … 124 125 125 126 ! Dynamic variables 126 real :: vcov(ip1jm,llm), ucov(ip1jmp1,llm) ! vents covariants 127 real :: teta(ip1jmp1,llm) ! temperature potentielle 128 real, dimension(:,:,:), allocatable :: q ! champs advectes 129 real :: ps(ip1jmp1) ! pression au sol 130 real, dimension(:), allocatable :: ps_start_GCM !(ngrid) pression au sol 131 real, dimension(:,:), allocatable :: ps_timeseries !(ngrid x timelen) ! pression au sol instantannées 132 real :: masse(ip1jmp1,llm) ! masse d'air 133 real :: phis(ip1jmp1) ! geopotentiel au sol 127 real, dimension(ip1jm,llm) :: vcov ! vents covariants 128 real, dimension(ip1jmp1,llm) :: ucov ! vents covariants 129 real, dimension(ip1jmp1,llm) :: teta ! temperature potentielle 130 real, dimension(:,:,:), allocatable :: q ! champs advectes 131 real, dimension(ip1jmp1) :: ps ! pression au sol 132 real, dimension(:), allocatable :: ps_start_GCM ! (ngrid) pression au sol 133 real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées 134 real, dimension(ip1jmp1,llm) :: masse ! masse d'air 135 real, dimension(ip1jmp1) :: phis ! geopotentiel au sol 134 136 real :: time_0 135 137 136 138 ! Variables to read starfi.nc 137 character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" ! Name of the file used for initialsing the PEM139 character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" ! Name of the file used for initialsing the PEM 138 140 character*2 str2 139 integer :: ncid, varid, status! Variable for handling opening of files141 integer :: ncid, varid, status ! Variable for handling opening of files 140 142 integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files 141 143 integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files … … 143 145 144 146 ! Variables to read starfi.nc and write restartfi.nc 145 real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAME and written in restartfi 146 real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi 147 real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart 148 real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart 149 real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi 147 real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAME and written in restartfi 148 real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi 149 real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi 150 150 real :: Total_surface ! Total surface of the planet 151 151 152 152 ! Variables for h2o_ice evolution 153 real :: ini_surf_h2o 154 real :: ini_surf_co2 155 real :: global_ave_press_GCM 156 real :: global_ave_press_old 157 real :: global_ave_press_new 158 real, dimension(:,:), allocatable :: zplev_new 159 real, dimension(:,:), allocatable :: zplev_gcm 160 real, dimension(:,:,:), allocatable :: zplev_new_timeseries 161 real, dimension(:,:,:), allocatable :: zplev_old_timeseries 162 logical :: STOPPING_water ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached?163 logical :: STOPPING_1_water ! Logical: is there still water ice to sublimate?164 logical :: STOPPING_co2 ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached?165 logical :: STOPPING_pressure ! Logical: is the criterion (% of change in the surface pressure) reached?166 integer :: criterion_stop 167 real, save :: A , B, mmean 168 real, allocatable :: vmr_co2_gcm(:,:)! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3]169 real, allocatable :: vmr_co2_pem_phys(:,:)! Physics x Times co2 volume mixing ratio used in the PEM170 real, allocatable :: q_co2_PEM_phys(:,:)! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg]171 real, allocatable :: q_h2o_PEM_phys(:,:)! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]172 integer :: timelen 173 real :: ave 174 real, allocatable :: p(:,:)! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)175 real :: extra_mass 153 real :: ini_surf_h2o ! Initial surface of sublimating h2o ice 154 real :: ini_surf_co2 ! Initial surface of sublimating co2 ice 155 real :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa] 156 real :: global_ave_press_old ! constant: Global average pressure of initial/previous time step 157 real :: global_ave_press_new ! constant: Global average pressure of current time step 158 real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2] 159 real, dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2] 160 real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] 161 real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step 162 logical :: STOPPING_water ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached? 163 logical :: STOPPING_1_water ! Logical: is there still water ice to sublimate? 164 logical :: STOPPING_co2 ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached? 165 logical :: STOPPING_pressure ! Logical: is the criterion (% of change in the surface pressure) reached? 166 integer :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param 167 real, save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] 168 real, dimension(:,:), allocatable :: vmr_co2_gcm ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3] 169 real, dimension(:,:), allocatable :: vmr_co2_pem_phys ! Physics x Times co2 volume mixing ratio used in the PEM 170 real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg] 171 real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg] 172 integer :: timelen ! # time samples 173 real :: ave ! intermediate varibale to compute average 174 real, dimension(:,:), allocatable :: p ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) 175 real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 176 176 177 177 ! Variables for slopes 178 real, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field 179 real, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field 180 real, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field 181 real, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field 182 real, dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field 183 real, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field 184 real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field 185 real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field 186 real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point x slope field: Tendency of evolution of perenial co2 ice over a year178 real, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field: minimum of co2 ice at each point for the first year [kg/m^2] 179 real, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field: minimum of co2 ice at each point for the second year [kg/m^2] 180 real, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field: minimum of water ice at each point for the first year [kg/m^2] 181 real, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field: minimum of water ice at each point for the second year [kg/m^2] 182 real, dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field: co2 ice given by the GCM [kg/m^2] 183 real, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field: Logical array indicating sublimating point of co2 ice 184 real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field: Logical array indicating if there is water ice at initial state 185 real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field: Logical array indicating if there is co2 ice at initial state 186 real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point x slope field: Tendency of evolution of perenial co2 ice over a year 187 187 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 188 real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical point x slope field: Tendency of evolution of perenial h2o ice189 real, dimension(:,:), allocatable :: flag_co2flow (:,:)! (ngrid,nslope): Flag where there is a CO2 glacier flow190 real, dimension(:), allocatable :: flag_co2flow_mesh (:)! (ngrid) : Flag where there is a CO2 glacier flow191 real, dimension(:,:), allocatable :: flag_h2oflow (:,:)! (ngrid,nslope): Flag where there is a H2O glacier flow192 real, dimension(:), allocatable :: flag_h2oflow_mesh (:)! (ngrid) : Flag where there is a H2O glacier flow188 real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical point x slope field: Tendency of evolution of perenial h2o ice 189 real, dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow 190 real, dimension(:), allocatable :: flag_co2flow_mesh ! (ngrid) : Flag where there is a CO2 glacier flow 191 real, dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow 192 real, dimension(:), allocatable :: flag_h2oflow_mesh ! (ngrid) : Flag where there is a H2O glacier flow 193 193 194 194 ! Variables for surface and soil 195 real, allocatable :: tsurf_ave(:,:) ! Physic x SLOPE field : Averaged Surface Temperature [K] 196 real, allocatable :: tsoil_ave(:,:,:) ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K] 197 real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K] 198 real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] 199 real, allocatable :: tsoil_GCM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] 200 real, allocatable :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] 201 real, allocatable :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] 202 real, allocatable :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K] 203 real, allocatable :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] 204 real, allocatable :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] 205 real, allocatable :: watersurf_density_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3] 206 real, allocatable :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] 207 real, allocatable :: watersoil_density_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3] 208 real, allocatable :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K] 209 real, allocatable :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] 210 real, allocatable :: delta_h2o_adsorbded(:) ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] 211 real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] 212 real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] 213 logical :: bool_sublim ! logical to check if there is sublimation or not 214 real,allocatable :: porefillingice_thickness_prev_iter(:,:) ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] 215 real,allocatable :: delta_h2o_icetablesublim(:) ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] 195 real, dimension(:,:), allocatable :: tsurf_ave ! Physic x SLOPE field : Averaged Surface Temperature [K] 196 real, dimension(:,:,:), allocatable :: tsoil_ave ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K] 197 real, dimension(:,:,:), allocatable :: tsurf_GCM_timeseries ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K] 198 real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] 199 real, dimension(:,:,:,:), allocatable :: tsoil_GCM_timeseries ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] 200 real, dimension(:,:), allocatable :: tsurf_ave_yr1 ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] 201 real, dimension(:,:), allocatable :: TI_locslope ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] 202 real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: intermediate when computing Tsoil [K] 203 real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] 204 real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] 205 real, dimension(:,:), allocatable :: watersurf_density_ave ! Physic x Slope, water surface density, yearly averaged [kg/m^3] 206 real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] 207 real, dimension(:,:,:), allocatable :: watersoil_density_PEM_ave ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3] 208 real, dimension(:,:), allocatable :: Tsurfave_before_saved ! Surface temperature saved from previous time step [K] 209 real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] 210 real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] 211 real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] 212 real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] 213 logical :: bool_sublim ! logical to check if there is sublimation or not 214 real, dimension(:,:), allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] 215 real, dimension(:), allocatable :: delta_h2o_icetablesublim(:) ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] 216 216 217 ! Some variables for the PEM run 217 218 real, parameter :: year_step = 1 ! timestep for the pem … … 221 222 integer :: n_myear ! Maximum number of Martian years of the chained simulations 222 223 real :: timestep ! timestep [s] 223 real :: watercap_sum ! total mass of water cap [kg/m^2] 224 real :: water_sum ! total mass of water in the mesh [kg/m^2] 224 real :: watercap_sum ! total mass of water cap [kg/m^2] 225 real :: water_sum ! total mass of water in the mesh [kg/m^2] 225 226 226 227 #ifdef CPP_STD 227 real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice228 real :: albedo_h2o_frost ! albedo of h2o frost229 real, allocatable :: tsurf_read_generic(:)! Temporary variable to do the subslope transfert dimensiion when reading form generic230 real, allocatable :: qsurf_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic231 real, allocatable :: tsoil_read_generic(:,:)! Temporary variable to do the subslope transfert dimensiion when reading form generic232 real, allocatable :: emis_read_generic(:)! Temporary variable to do the subslope transfert dimensiion when reading form generic233 real, allocatable :: albedo_read_generic(:,:)! Temporary variable to do the subslope transfert dimensiion when reading form generic234 real, allocatable :: tsurf(:,:)! Subslope variable, only needed in the GENERIC case235 real, allocatable :: qsurf(:,:,:)! Subslope variable, only needed in the GENERIC case236 real, allocatable :: tsoil(:,:,:)! Subslope variable, only needed in the GENERIC case237 real, allocatable :: emis(:,:)! Subslope variable, only needed in the GENERIC case238 real, allocatable :: watercap(:,:)! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model239 logical, allocatable :: WATERCAPTAG(:)! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model240 real, allocatable :: albedo(:,:,:)! Subslope variable, only needed in the GENERIC case241 real, allocatable :: inertiesoil(:,:,:)! Subslope variable, only needed in the GENERIC case228 real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice 229 real :: albedo_h2o_frost ! albedo of h2o frost 230 real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic 231 real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic 232 real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic 233 real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic 234 real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic 235 real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case 236 real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case 237 real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case 238 real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case 239 real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model 240 logical, dimension(:), allocatable :: WATERCAPTAG ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model 241 real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case 242 real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case 242 243 #endif 243 244 244 245 #ifdef CPP_1D 245 integer :: nsplit_phys 246 integer, parameter :: jjm_value = jjm - 1 247 integer :: ierr 246 integer :: nsplit_phys 247 integer, parameter :: jjm_value = jjm - 1 248 integer :: ierr 249 250 ! Dummy variables to use the subroutine 'init_testphys1d' 251 real, dimension(:,:), allocatable :: q_tmp ! Temporary tracer mixing ratio (e.g. kg/kg) 252 logical :: startfiles_1D, therestart1D, therestartfi 253 integer :: ndt, day0 254 real :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau 255 real, dimension(:), allocatable :: zqsat, qsurf_tmp 256 real, dimension(:,:), allocatable :: dq, dqdyn 257 real, dimension(nlayer) :: play, w 258 real, dimension(nlayer + 1) :: plev, q2_tmp 259 real, dimension(nsoilmx) :: tsoil_tmp 260 real, dimension(1) :: tsurf_tmp, emis_tmp 261 real, dimension(1,1) :: albedo_tmp 248 262 #else 249 integer, parameter :: jjm_value = jjm 263 integer, parameter :: jjm_value = jjm 264 real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart 265 real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart 250 266 #endif 251 267 … … 262 278 #endif 263 279 280 ! Some constants 264 281 day_ini = 0 ! test 265 282 time_phys = 0. ! test 266 267 ! Some constants268 283 ngrid = ngridmx 269 nlayer = llm270 271 284 A = (1/m_co2 - 1/m_noco2) 272 285 B = 1/m_noco2 273 274 286 year_day = 669 275 287 daysec = 88775. … … 287 299 nq = nqtot 288 300 allocate(q(ip1jmp1,llm,nqtot)) 301 allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid)) 289 302 #else 290 ! load tracer names from file 'traceur.def' 291 open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr) 292 if (ierr /= 0) then 293 write(*,*) 'Cannot find required file "traceur.def"' 294 write(*,*) ' If you want to run with tracers, I need it' 295 write(*,*) ' ... might as well stop here ...' 296 stop 297 else 298 write(*,*) "pem1d: Reading file traceur.def" 299 ! read number of tracers: 300 read(90,*,iostat = ierr) nq 301 write(*,*) "nq",nq 302 nqtot = nq ! set value of nqtot (in infotrac module) as nq 303 if (ierr /= 0) then 304 write(*,*) "pem1d: error reading number of tracers" 305 write(*,*) " (first line of traceur.def) " 306 stop 307 endif 308 if (nq < 1) then 309 write(*,*) "pem1d: error number of tracers" 310 write(*,*) "is nq=",nq," but must be >=1!" 311 stop 312 endif 313 endif 314 nq = nqtot 303 allocate(longitude(1),latitude(1),cell_area(1)) 304 call init_testphys1d(.true.,ngrid,nlayer,610.,nq,q_tmp,time_0,ps(1),ucov,vcov,teta,startfiles_1D,therestart1D,therestartfi, & 305 ndt,ptif,pks,dtphys,zqsat,qsurf_tmp,dq,dqdyn,day0,day,tsurf_tmp,gru,grv,w,q2_tmp,play,plev,tsoil_tmp, & 306 albedo_tmp,emis_tmp,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau) 307 ps(2) = ps(1) 315 308 allocate(q(ip1jmp1,llm,nqtot)) 316 allocate(ap(nlayer + 1)) 317 allocate(bp(nlayer + 1)) 318 call init_pem1d(llm,nqtot,ucov,vcov,teta,q,ps,time_0,ap,bp) 319 pi = 2.*asin(1.) 320 g = 3.72 309 q(ip1jmp1,:,:) = q_tmp(:,:) 310 deallocate(q_tmp,qsurf_tmp) 321 311 nsplit_phys = 1 322 312 #endif … … 356 346 ! Here we simply read them in the startfi_evol.nc file 357 347 status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) 358 359 allocate(longitude(ngrid))360 allocate(latitude(ngrid))361 allocate(cell_area(ngrid))362 348 363 349 status = nf90_inq_varid(ncid,"longitude",lonvarid) … … 414 400 qsurf_read_generic,cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, & 415 401 tslab, tsea_ice,sea_ice) 416 call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 402 call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 417 403 418 404 nslope = 1 … … 452 438 endif 453 439 if (noms(nnq) == "co2") igcm_co2 = nnq 454 enddo 440 enddo 455 441 r = 8.314511*1000./mugaz 456 442 … … 473 459 allocate(flag_h2oflow_mesh(ngrid)) 474 460 475 flag_co2flow(:,:) = 0 461 flag_co2flow(:,:) = 0 476 462 flag_co2flow_mesh(:) = 0 477 flag_h2oflow(:,:) = 0 463 flag_h2oflow(:,:) = 0 478 464 flag_h2oflow_mesh(:) = 0 479 465 … … 614 600 global_ave_press_old = 0. 615 601 do i = 1,ngrid 616 602 global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface 617 603 enddo 618 604 … … 656 642 write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded 657 643 endif ! adsorption 658 deallocate(tsurf_ave_yr1) 644 deallocate(tsurf_ave_yr1) 659 645 660 646 !------------------------ … … 678 664 !------------------------ 679 665 ! II Run 680 ! II_a Update pressure, ice and tracers 666 ! II_a Update pressure, ice and tracers 681 667 !------------------------ 682 668 year_iter = 0 … … 688 674 do i = 1,ngrid 689 675 do islope = 1,nslope 690 global_ave_press_new = global_ave_press_new -g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface691 enddo 692 enddo 693 call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)694 676 global_ave_press_new = global_ave_press_new - g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface 677 enddo 678 enddo 679 call writediagfi(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new) 680 695 681 if (adsorption_pem) then 696 682 do i = 1,ngrid … … 739 725 740 726 do ig = 1,ngrid 741 do t = 1, 727 do t = 1,timelen 742 728 q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ & 743 729 (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & … … 767 753 call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope) 768 754 769 do islope =1,nslope755 do islope = 1,nslope 770 756 write(str2(1:2),'(i2.2)') islope 771 call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))772 call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))773 call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))774 call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))757 call writediagfi(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope)) 758 call writediagfi(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope)) 759 call writediagfi(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) 760 call writediagfi(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) 775 761 enddo 776 762 … … 783 769 if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, & 784 770 global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh) 785 771 786 772 write(*,*) "H2O glacier flows" 787 773 … … 790 776 do islope = 1,nslope 791 777 write(str2(1:2),'(i2.2)') islope 792 call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))793 call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))794 call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))778 call writediagfi(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) 779 call writediagfi(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) 780 call writediagfi(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) 795 781 enddo 796 782 … … 835 821 else 836 822 albedo(ig,1,islope) = albedodat(ig) 837 albedo(ig,2,islope) = albedodat(ig) 823 albedo(ig,2,islope) = albedodat(ig) 838 824 emis(ig,islope) = emissiv 839 825 endif … … 862 848 863 849 do t = 1,timelen 864 tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:)) 850 tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:)) 865 851 enddo 866 852 ! for the start … … 873 859 do islope = 1,nslope 874 860 write(str2(1:2),'(i2.2)') islope 875 call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))861 call writediagfi(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) 876 862 enddo 877 863 … … 896 882 do isoil = 1,nsoilmx_PEM 897 883 watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r) 898 if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1) 884 if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1) 899 885 enddo 900 886 enddo … … 1001 987 endif 1002 988 1003 global_ave_press_old =global_ave_press_new989 global_ave_press_old = global_ave_press_new 1004 990 1005 991 enddo ! big time iteration loop of the pem … … 1037 1023 water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 1038 1024 enddo 1039 if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. 1025 if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. 1040 1026 if (water_sum > threshold_water_frost2perenial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1 1041 1027 watercaptag(ig) = .true. … … 1071 1057 if (soil_pem) then 1072 1058 call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) 1073 tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen) 1059 tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen) 1074 1060 endif 1075 1061 1076 1062 ! III_a.4 Pressure (for start) 1077 do i = 1,ip1jmp1 1078 ps(i) = ps(i)*global_ave_press_new/global_ave_press_GCM 1079 enddo 1080 1081 do i = 1,ngrid 1082 ps_start_GCM(i) = ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM 1083 enddo 1063 ps(:) = ps(:)*global_ave_press_new/global_ave_press_GCM 1064 ps_start_GCM(:) = ps_start_GCM(:)*global_ave_press_new/global_ave_press_GCM 1084 1065 1085 1066 ! III_a.5 Tracer (for start) … … 1087 1068 1088 1069 do l = 1,nlayer + 1 1089 do ig = 1,ngrid 1090 zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) 1091 enddo 1070 zplev_new(:,l) = ap(l) + bp(l)*ps_start_GCM(:) 1092 1071 enddo 1093 1072 … … 1104 1083 do ig = 1,ngrid 1105 1084 q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) & 1106 + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/ &1085 + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/ & 1107 1086 (zplev_new(ig,l) - zplev_new(ig,l + 1)) 1108 1087 enddo … … 1117 1096 do l = 1,llm - 1 1118 1097 if (q(ig,l,nnq) > 1 .and. (noms(nnq) /= "dust_number") .and. (noms(nnq) /= "ccn_number") .and. (noms(nnq) /= "stormdust_number") .and. (noms(nnq) /= "topdust_number")) then 1119 extra_mass =(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))1120 q(ig,l,nnq) =1.1121 q(ig,l +1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))1098 extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1)) 1099 q(ig,l,nnq) = 1. 1100 q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2)) 1122 1101 write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number" 1123 1102 endif … … 1146 1125 write(*,*) "restart_evol.nc has been written" 1147 1126 #else 1148 do nnq = 1, nqtot 1149 call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf) 1150 enddo 1127 call writerestart1D('restart1D_evol.txt',ps(1),tsurf,nlayer,teta,ucov,vcov,nq,noms,qsurf,q) 1128 write(*,*) "restart1D_evol.txt has been written" 1151 1129 #endif 1152 1130 … … 1198 1176 1199 1177 END PROGRAM pem 1200 -
trunk/LMDZ.MARS/deftank/pem/launch_pem.sh
r3050 r3065 80 80 if [ ! -d "diags" ]; then 81 81 mkdir diags 82 fi83 if [ ! -d "profiles" ]; then84 mkdir profiles85 82 fi 86 83 … … 98 95 ((iGCM = ($iPEM - 1)*$nGCM + 1)) 99 96 cp startfi.nc starts/ 100 for file in profile_*; do 101 cp $file profiles/${file}_0102 done 97 if [ -f "start.nc" ]; then 98 cp start.nc starts/ 99 fi 103 100 104 101 # Create a temporary file to manage years of the chained simulation … … 138 135 cp restart.nc starts/restart${iGCM}.nc 139 136 mv restart.nc start.nc 140 fi 141 if [ -f "restart1D.txt" ]; then 137 elif [ -f "restart1D.txt" ]; then 142 138 cp restart1D.txt starts/restart1D${iGCM}.txt 143 139 mv restart1D.txt start1D.txt 144 fi145 if ls profile_out_* 1> /dev/null 2>&1; then146 for file in profile_out_*; do147 cp $file profiles/${file/_out/}_${iGCM}148 mv $file ${file/_out/}149 done150 140 fi 151 141 ((iGCM++)) … … 171 161 if [ -f "start.nc" ]; then 172 162 cp start.nc start_evol.nc 173 fi 174 #if [ -f "start1D.txt" ]; then 175 # cp start1D.txt start1D.txt 176 #fi 163 elif [ -f "start1D.txt" ]; then 164 cp start1D.txt start1D_evol.txt 165 fi 177 166 ./$exePEM > out_runPEM${iPEM} 2>&1 178 167 if [ ! -f "restartfi_evol.nc" ]; then # Check if run ended abnormally … … 183 172 mv out_runPEM${iPEM} out_PEM/run${iPEM} 184 173 mv diagfi.nc diags/diagfi_PEM${iPEM}.nc 174 cp restartfi_PEM.nc starts/startfi_PEM${iPEM}.nc 175 mv restartfi_PEM.nc startfi_PEM.nc 185 176 cp restartfi_evol.nc starts/startfi_postPEM${iPEM}.nc 186 177 mv restartfi_evol.nc startfi.nc 187 178 if [ -f "restart_evol.nc" ]; then 188 cp restart_evol.nc starts/restart ${iGCM}.nc179 cp restart_evol.nc starts/restart_postPEM${iPEM}.nc 189 180 mv restart_evol.nc start.nc 190 fi 191 cp restartfi_PEM.nc starts/startfi_PEM${iPEM}.nc 192 mv restartfi_PEM.nc startfi_PEM.nc 193 if ls profile_out_* 1> /dev/null 2>&1; then 194 for file in profile_out_*; do 195 cp $file profiles/${file/_out/}PEM_${iPEM} 196 mv $file ${file/_out/} 197 done 181 elif [ -f "restart1D_evol.txt" ]; then 182 cp restart1D_evol.txt starts/restart1D_postPEM${iPEM}.txt 183 mv restart1D_evol.txt start1D.txt 198 184 fi 199 185 ((iPEM++)) -
trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
r3060 r3065 5 5 contains 6 6 7 SUBROUTINE init_testphys1d( startfiles_1D,therestart1D,therestartfi,ngrid,nlayer,odpref,nq,ndt,ptif,pks,dttestphys, &8 q,zqsat,qsurf,dq,dqdyn,day0,day,time,psurf,tsurf,gru,grv,u,v,w,q2,play,plev,tsoil,temp,&7 SUBROUTINE init_testphys1d(pem1d,ngrid,nlayer,odpref,nq,q,time,psurf,u,v,temp,startfiles_1D,therestart1D,therestartfi, & 8 ndt,ptif,pks,dttestphys,zqsat,qsurf,dq,dqdyn,day0,day,tsurf,gru,grv,w,q2,play,plev,tsoil, & 9 9 albedo,emis,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau) 10 10 … … 47 47 ! Arguments 48 48 !======================================================================= 49 logical, intent(in) :: startfiles_1D, therestart1D, therestartfi ! Use of "start1D.txt" and "startfi.nc" files 50 integer, intent(in) :: ngrid, nlayer 51 real, intent(in) :: odpref ! DOD reference pressure (Pa) 49 integer, intent(in) :: ngrid, nlayer 50 real, intent(in) :: odpref ! DOD reference pressure (Pa) 51 logical, intent(in) :: pem1d ! If initialization for the 1D PEM 52 52 53 53 integer, intent(inout) :: nq 54 54 55 real, dimension(:,:), allocatable, intent(out) :: q ! tracer mixing ratio (e.g. kg/kg) 56 real, intent(out) :: time ! time (0<time<1; time=0.5 at noon) 57 real, intent(out) :: psurf ! Surface pressure 58 real, dimension(nlayer), intent(out) :: u, v ! zonal, meridional wind 59 real, dimension(nlayer), intent(out) :: temp ! temperature at the middle of the layers 60 logical, intent(out) :: startfiles_1D, therestart1D, therestartfi ! Use of starting files for 1D 55 61 integer, intent(out) :: ndt 56 62 real, intent(out) :: ptif, pks 57 real, intent(out) :: dttestphys ! testphys1d timestep 58 real, dimension(:,:), allocatable, intent(out) :: q ! tracer mixing ratio (e.g. kg/kg) 59 real, dimension(:), allocatable, intent(out) :: zqsat ! useful for (atm_wat_profile=2) 60 real, dimension(:), allocatable, intent(out) :: qsurf ! tracer surface budget (e.g. kg.m-2) 61 real, dimension(:,:), allocatable, intent(out) :: dq, dqdyn ! Physical and dynamical tandencies 62 integer, intent(out) :: day0 ! initial (sol ; =0 at Ls=0) and final date 63 real, intent(out) :: day ! date during the run 64 real, intent(out) :: time ! time (0<time<1; time=0.5 at noon) 65 real, intent(out) :: psurf ! Surface pressure 66 real, dimension(1), intent(out) :: tsurf ! Surface temperature 67 real, intent(out) :: gru, grv ! prescribed "geostrophic" background wind 68 real, dimension(nlayer), intent(out) :: u, v, w ! zonal, meridional wind 69 real, dimension(nlayer + 1), intent(out) :: q2 ! Turbulent Kinetic Energy 70 real, dimension(nlayer), intent(out) :: play ! Pressure at the middle of the layers (Pa) 71 real, dimension(nlayer + 1), intent(out) :: plev ! intermediate pressure levels (pa) 72 real, dimension(nsoilmx), intent(out) :: tsoil ! subsurface soil temperature (K) 73 real, dimension(nlayer), intent(out) :: temp ! temperature at the middle of the layers 74 real, dimension(1,1), intent(out) :: albedo ! surface albedo 75 real, dimension(1), intent(out) :: emis ! surface layer 63 real, intent(out) :: dttestphys ! testphys1d timestep 64 real, dimension(:), allocatable, intent(out) :: zqsat ! useful for (atm_wat_profile=2) 65 real, dimension(:), allocatable, intent(out) :: qsurf ! tracer surface budget (e.g. kg.m-2) 66 real, dimension(:,:), allocatable, intent(out) :: dq, dqdyn ! Physical and dynamical tandencies 67 integer, intent(out) :: day0 ! initial (sol ; =0 at Ls=0) and final date 68 real, intent(out) :: day ! date during the run 69 real, dimension(1), intent(out) :: tsurf ! Surface temperature 70 real, intent(out) :: gru, grv ! prescribed "geostrophic" background wind 71 real, dimension(nlayer), intent(out) :: w ! "Dummy wind" in 1D 72 real, dimension(nlayer + 1), intent(out) :: q2 ! Turbulent Kinetic Energy 73 real, dimension(nlayer), intent(out) :: play ! Pressure at the middle of the layers (Pa) 74 real, dimension(nlayer + 1), intent(out) :: plev ! intermediate pressure levels (pa) 75 real, dimension(nsoilmx), intent(out) :: tsoil ! subsurface soil temperature (K) 76 real, dimension(1,1), intent(out) :: albedo ! surface albedo 77 real, dimension(1), intent(out) :: emis ! surface layer 76 78 real, dimension(1), intent(out) :: latitude, longitude, cell_area 77 real, intent(out) :: atm_wat_profile, atm_wat_tau 79 real, intent(out) :: atm_wat_profile, atm_wat_tau ! Force atmospheric water profiles 78 80 79 81 !======================================================================= … … 89 91 real, dimension(:,:), allocatable :: psdyn 90 92 91 ! RV & JBC: Use of "start1D.txt" and "startfi.nc" files93 ! RV & JBC: Use of starting files for 1D 92 94 logical :: found 93 95 character(len = 30) :: header … … 102 104 integer :: ifils, ipere, generation, ierr0 103 105 character(len = 30), dimension(:), allocatable :: tnom_transp ! transporting fluid short name 104 character(len = 80) :: line ! to store a line of text106 character(len = 80) :: line ! to store a line of text 105 107 logical :: continu, there 106 108 … … 113 115 real :: flux_geo_tmp 114 116 117 ! JBC: To initialize the 1D PEM 118 character(:), allocatable :: start1Dname, startfiname ! Name of starting files for 1D 119 115 120 !======================================================================= 116 121 ! Code 117 122 !======================================================================= 123 if (.not. pem1d) then 124 start1Dname = 'start1D.txt' 125 startfiname = 'startfi.nc' 126 startfiles_1D = .false. 127 !------------------------------------------------------ 128 ! Loading run parameters from "run.def" file 129 !------------------------------------------------------ 130 ! check if 'run.def' file is around. Otherwise reading parameters 131 ! from callphys.def via getin() routine won't work. 132 inquire(file = 'run.def',exist = there) 133 if (.not. there) then 134 write(*,*) 'Cannot find required file "run.def"' 135 write(*,*) ' (which should contain some input parameters along with the following line: INCLUDEDEF=callphys.def)' 136 write(*,*) ' ... might as well stop here ...' 137 stop 138 endif 139 140 write(*,*)'Do you want to use starting files?' 141 call getin("startfiles_1D",startfiles_1D) 142 write(*,*) " startfiles_1D = ", startfiles_1D 143 else 144 start1dname = 'start1D_evol.txt' 145 startfiname = 'startfi_evol.nc' 146 startfiles_1D = .true. 147 endif 148 149 therestart1D = .false. 150 therestartfi = .false. 151 inquire(file = start1Dname,exist = therestart1D) 152 if (startfiles_1D .and. .not. therestart1D) then 153 write(*,*) 'There is no "'//start1Dname//'" file!' 154 if (.not. pem1d) then 155 write(*,*) 'Initialization is done with default values.' 156 else 157 write(*,*) 'Initialization cannot be done for the 1D PEM.' 158 stop 159 endif 160 endif 161 inquire(file = startfiname,exist = therestartfi) 162 if (.not. therestartfi) then 163 write(*,*) 'There is no "'//startfiname//'" file!' 164 if (.not. pem1d) then 165 write(*,*) 'Initialization is done with default values.' 166 else 167 write(*,*) 'Initialization cannot be done for the 1D PEM.' 168 stop 169 endif 170 endif 118 171 119 172 !------------------------------------------------------ … … 186 239 endif 187 240 endif 241 188 242 ! allocate arrays: 189 243 allocate(tname(nq),q(nlayer,nq),zqsat(nlayer),qsurf(nq)) … … 245 299 else 246 300 do iq = 1,nq 247 open(3,file = 'start1D.txt',status = "old",action = "read")301 open(3,file = start1Dname,status = "old",action = "read") 248 302 read(3,*) header, qsurf(iq),(q(ilayer,iq), ilayer = 1,nlayer) 249 303 if (trim(tname(iq)) /= trim(header)) then 250 write(*,*) 'Tracer names not compatible for initialization with " start1D.txt"!'304 write(*,*) 'Tracer names not compatible for initialization with "'//trim(start1Dname)//'"!' 251 305 stop 252 306 endif … … 254 308 endif 255 309 256 call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,COMM_LMDZ) 310 #ifdef CPP_XIOS 311 call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,COMM_LMDZ) 312 #else 313 call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,1) 314 #endif 257 315 258 316 ! Date and local time at beginning of run … … 261 319 ! Date (in sols since spring solstice) at beginning of run 262 320 day0 = 0 ! default value for day0 263 write(*,*) 'Initial date (in martian sols 321 write(*,*) 'Initial date (in martian sols; =0 at Ls=0)?' 264 322 call getin("day0",day0) 265 323 day = float(day0) … … 272 330 time = time/24. ! convert time (hours) to fraction of sol 273 331 else 274 call open_startphy( "startfi.nc")332 call open_startphy(startfiname) 275 333 call get_var("controle",tab_cntrl,found) 276 334 if (.not. found) then … … 545 603 zlay(:) = -200.*r*log(play(:)/plev(1))/g 546 604 547 548 605 ! Initialize temperature profile 549 606 ! ------------------------------ -
trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90
r3060 r3065 1 1 PROGRAM testphys1d 2 2 3 use ioipsl_getincom, only: getin ! To use 'getin'4 3 use comsoil_h, only: inertiedat, inertiesoil, nsoilmx 5 4 use surfdat_h, only: albedodat, perenial_co2ice, watercap … … 100 99 character(len = 44) :: txt 101 100 102 ! RV & JBC: Use of "start1D.txt" and "startfi.nc" files101 ! RV & JBC: Use of starting files for 1D 103 102 logical :: startfiles_1D, therestart1D, therestartfi 104 103 … … 121 120 !call initcomgeomphy 122 121 123 !------------------------------------------------------ 124 ! Loading run parameters from "run.def" file 125 !------------------------------------------------------ 126 ! check if 'run.def' file is around. Otherwise reading parameters 127 ! from callphys.def via getin() routine won't work. 128 inquire(file = 'run.def',exist = there) 129 if (.not. there) then 130 write(*,*) 'Cannot find required file "run.def"' 131 write(*,*) ' (which should contain some input parameters along with the following line: INCLUDEDEF=callphys.def)' 132 write(*,*) ' ... might as well stop here ...' 133 stop 134 endif 135 136 write(*,*)'Do you want to use "start1D.txt" and "startfi.nc" files?' 137 startfiles_1D = .false. 138 therestart1D = .false. 139 therestartfi = .false. 140 call getin("startfiles_1D",startfiles_1D) 141 write(*,*) " startfiles_1D = ", startfiles_1D 142 143 if (startfiles_1D) then 144 inquire(file = 'start1D.txt',exist = therestart1D) 145 if (.not. therestart1D) then 146 write(*,*) 'There is no "start1D.txt" file!' 147 write(*,*) 'Initialization is done with default values.' 148 endif 149 inquire(file = 'startfi.nc',exist = therestartfi) 150 if (.not. therestartfi) then 151 write(*,*) 'There is no "startfi.nc" file!' 152 write(*,*) 'Initialization is done with default values.' 153 endif 154 endif 155 156 call init_testphys1d(startfiles_1D,therestart1D,therestartfi,ngrid,nlayer,odpref,nq,ndt,ptif,pks,dttestphys, & 157 q,zqsat,qsurf,dq,dqdyn,day0,day,time,psurf,tsurf,gru,grv,u,v,w,q2,play,plev,tsoil,temp, & 122 call init_testphys1d(.false.,ngrid,nlayer,odpref,nq,q,time,psurf,u,v,temp,startfiles_1D,therestart1D,therestartfi, & 123 ndt,ptif,pks,dttestphys,zqsat,qsurf,dq,dqdyn,day0,day,tsurf,gru,grv,w,q2,play,plev,tsoil, & 158 124 albedo,emis,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau) 159 125 … … 265 231 266 232 ! Writing the "restart1D.txt" file for the next run 267 if (startfiles_1D) call writerestart1D( psurf,tsurf,nlayer,temp,u,v,nq,noms,qsurf,q)233 if (startfiles_1D) call writerestart1D('restart1D.txt',psurf,tsurf,nlayer,temp,u,v,nq,noms,qsurf,q) 268 234 269 235 write(*,*) "testphys1d: everything is cool." -
trunk/LMDZ.MARS/libf/phymars/dyn1d/writerestart1D.F90
r3051 r3065 1 SUBROUTINE writerestart1D( psurf,tsurf,nlayer,temp,u,v,nq,qnames,qsurf,q)1 SUBROUTINE writerestart1D(filename,psurf,tsurf,nlayer,temp,u,v,nq,qnames,qsurf,q) 2 2 3 3 implicit none 4 4 5 5 ! Arguments 6 integer, intent(in) :: nlayer, nq 7 real, intent(in) :: psurf, tsurf 8 real, dimension(nlayer), intent(in) :: temp, u, v 9 real, dimension(nlayer,nq), intent(in) :: q 10 real, dimension(nq), intent(in) :: qsurf 6 character(len = *), intent(in) :: filename 7 integer, intent(in) :: nlayer, nq 8 real, intent(in) :: psurf, tsurf 9 real, dimension(nlayer), intent(in) :: temp, u, v 10 real, dimension(nlayer,nq), intent(in) :: q 11 real, dimension(nq), intent(in) :: qsurf 11 12 character(len = *), dimension(nq), intent(in) :: qnames 12 13 … … 15 16 16 17 ! Write the data needed for a restart in "restart1D.txt" 17 open(1,file = 'restart1D.txt',status = "replace",action = "write")18 open(1,file = filename,status = "replace",action = "write") 18 19 do iq = 1,nq 19 20 write(1,*) qnames(iq), qsurf(iq), (q(il,iq), il = 1,nlayer)
Note: See TracChangeset
for help on using the changeset viewer.
