Changeset 3065 for trunk/LMDZ.COMMON/libf
- Timestamp:
- Oct 2, 2023, 2:30:47 PM (22 months ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 1 deleted
- 2 edited
-
changelog.txt (modified) (1 diff)
-
init_pem1d_mod.F90 (deleted)
-
pem.F90 (modified) (31 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
Note: See TracChangeset
for help on using the changeset viewer.
