Changeset 3602 for trunk/LMDZ.COMMON/libf
- Timestamp:
- Jan 27, 2025, 5:44:22 PM (4 days ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 6 edited
-
changelog.txt (modified) (1 diff)
-
compute_soiltemp_mod.F90 (modified) (5 diffs)
-
ice_table_mod.F90 (modified) (5 diffs)
-
pem.F90 (modified) (1 diff)
-
pemetat0.F90 (modified) (9 diffs)
-
read_data_PCM_mod.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3599 r3602 557 557 == 23/01/2025 == JBC 558 558 Few small optimizations 559 560 == 27/01/2025 == JBC 561 - Complete rework of "read_data_PCM_mod.F90" to reduce as much as possible memory consumption which caused problems and the model failure. 562 - In case of very low surface shifting, checking if the index is inside the bounds of 'tsoil' and 'mlayer' and replace them by the value at surface if necessary. 563 - Few cleanings. -
trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90
r3553 r3602 243 243 ! ------ 244 244 integer :: ig, isoil, islope, ishift, ilag, iz 245 real :: a, z, zshift_surfloc 245 real :: a, z, zshift_surfloc, tsoil_minus, mlayer_minus 246 246 real, dimension(ngrid,nsoil,nslope) :: tsoil_old 247 247 … … 265 265 ! The "new soil" temperature is set to tsurf 266 266 tsoil(ig,:ishift,islope) = tsurf(ig,islope) 267 268 267 do isoil = ishift + 1,nsoil 269 268 ! Position in the old discretization of the depth … … 274 273 iz = iz + 1 275 274 enddo 275 if (iz == 0) then 276 tsoil_minus = tsurf(ig,islope) 277 mlayer_minus = 0. 278 else 279 tsoil_minus = tsoil_old(ig,iz,islope) 280 mlayer_minus = mlayer_PEM(iz - 1) 281 endif 276 282 ! Interpolation of the temperature profile from the old discretization 277 a = (tsoil_old(ig,iz + 1,islope) - tsoil_ old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))278 tsoil(ig,isoil,islope) = a*(z - mlayer_ PEM(iz - 1)) + tsoil_old(ig,iz,islope)283 a = (tsoil_old(ig,iz + 1,islope) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus) 284 tsoil(ig,isoil,islope) = a*(z - mlayer_minus) + tsoil_minus 279 285 enddo 280 286 endif … … 299 305 iz = iz + 1 300 306 enddo 307 if (iz == 0) then 308 tsoil_minus = tsurf(ig,islope) 309 mlayer_minus = 0. 310 else 311 tsoil_minus = tsoil_old(ig,iz,islope) 312 mlayer_minus = mlayer_PEM(iz - 1) 313 endif 301 314 ! The "new lag" temperature is set to the ice temperature just below 302 a = (tsoil_old(ig,iz + 1,islope) - tsoil_ old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))303 tsoil(ig,:ilag,islope) = a*(z - mlayer_ PEM(iz - 1)) + tsoil_old(ig,iz,islope)315 a = (tsoil_old(ig,iz + 1,islope) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus) 316 tsoil(ig,:ilag,islope) = a*(z - mlayer_minus) + tsoil_minus 304 317 endif 305 318 … … 318 331 iz = iz + 1 319 332 enddo 333 if (iz == 0) then 334 tsoil_minus = tsurf(ig,islope) 335 mlayer_minus = 0. 336 else 337 tsoil_minus = tsoil_old(ig,iz,islope) 338 mlayer_minus = mlayer_PEM(iz - 1) 339 endif 320 340 ! Interpolation of the temperature profile from the old discretization 321 a = (tsoil_old(ig,iz + 1,islope) - tsoil_ old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))322 tsoil(ig,isoil,islope) = a*(z - mlayer_ PEM(iz - 1)) + tsoil_old(ig,iz,islope)341 a = (tsoil_old(ig,iz + 1,islope) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus) 342 tsoil(ig,isoil,islope) = a*(z - mlayer_minus) + tsoil_minus 323 343 enddo 324 344 -
trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
r3571 r3602 46 46 47 47 !------------------------------------------------------------------------------------------------------ 48 SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_ beg,ice_table_thickness)48 SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_depth,ice_table_thickness) 49 49 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 50 50 !!! … … 72 72 ! Ouputs 73 73 ! ------ 74 real, dimension(ngrid,nslope), intent(out) :: ice_table_ beg! ice table depth [m]74 real, dimension(ngrid,nslope), intent(out) :: ice_table_depth ! ice table depth [m] 75 75 real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m] 76 76 ! Locals … … 84 84 ! Code 85 85 ! ---- 86 previous_icetable_depth = ice_table_ beg86 previous_icetable_depth = ice_table_depth 87 87 do ig = 1,ngrid 88 88 if (watercaptag(ig)) then 89 ice_table_ beg(ig,:) = 0.89 ice_table_depth(ig,:) = 0. 90 90 ice_table_thickness(ig,:) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) 91 91 else 92 92 do islope = 1,nslope 93 ice_table_ beg(ig,islope) = -1.93 ice_table_depth(ig,islope) = -1. 94 94 ice_table_thickness(ig,islope) = 0. 95 95 do isoil = 1,nsoil 96 96 diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope) 97 97 enddo 98 if (diff_rho(1) > 0 ) then ! ice is at the surface99 ice_table_ beg(ig,islope) = 0.98 if (diff_rho(1) > 0.) then ! ice is at the surface 99 ice_table_depth(ig,islope) = 0. 100 100 do isoilend = 2,nsoil - 1 101 101 if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then 102 call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend ),mlayer_PEM(isoilend + 1),ice_table_end)103 ice_table_thickness(ig,islope) = ice_table_end - ice_table_ beg(ig,islope)102 call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend - 1),mlayer_PEM(isoilend),ice_table_end) 103 ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope) 104 104 exit 105 105 endif … … 107 107 else 108 108 do isoil = 1,nsoil - 1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root 109 if (diff_rho(isoil) < 0 .and. diff_rho(isoil + 1) > 0.) then110 call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_ beg(ig,islope))109 if (diff_rho(isoil) < 0. .and. diff_rho(isoil + 1) > 0.) then 110 call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_depth(ig,islope)) 111 111 ! Now let's find the end of the ice table 112 112 ice_table_thickness(ig,islope) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) 113 113 do isoilend = isoil + 1,nsoil - 1 114 if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then115 call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend ),mlayer_PEM(isoilend + 1),ice_table_end)116 ice_table_thickness(ig,islope) = ice_table_end - ice_table_ beg(ig,islope)114 if (diff_rho(isoilend) > 0. .and. diff_rho(isoilend + 1) < 0.) then 115 call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend - 1),mlayer_PEM(isoilend),ice_table_end) 116 ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope) 117 117 exit 118 118 endif … … 128 128 do islope = 1,nslope 129 129 do ig = 1,ngrid 130 if (ice_table_ beg(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then130 if (ice_table_depth(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then 131 131 call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia) 132 132 stretch = (regolith_inertia(ig,islope)/wice_inertia)**2 133 ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_ beg(ig,islope) - &134 previous_icetable_depth(ig,islope) + (ice_table_ beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)135 ice_table_ beg(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch133 ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_depth(ig,islope) - & 134 previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch) 135 ice_table_depth(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch 136 136 endif 137 137 enddo -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3599 r3602 615 615 616 616 delta_h2o_icetablesublim = 0. 617 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, & 618 icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys, & 619 ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice,watersurf_density_avg, & 620 watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif) 617 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, & 618 tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice, & 619 watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif) 621 620 deallocate(tsurf_avg_yr1) 622 621 -
trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
r3571 r3602 7 7 !======================================================================= 8 8 9 SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness, & 10 ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global, & 11 watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys, & 12 d_h2oice,d_co2ice,co2_ice,h2o_ice,m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif) 9 SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, & 10 tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global,d_h2oice,d_co2ice,co2_ice,h2o_ice, & 11 watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif) 13 12 14 13 use iostart_PEM, only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length … … 210 209 enddo ! islope 211 210 212 write(*,*) 'PEMETAT0: T HERMAL INERTIAdone'211 write(*,*) 'PEMETAT0: TI done' 213 212 214 213 ! b. Special case for inertiedat, inertiedat_PEM … … 277 276 enddo 278 277 enddo ! islope 279 write(*,*) 'PEMETAT0: SOIL TEMPdone'278 write(*,*) 'PEMETAT0: TSOIL done' 280 279 281 280 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 367 366 !!! transition 368 367 delta = depth_breccia 369 TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia +1)-layer_PEM(index_breccia))/&370 (((delta -layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &371 ((layer_PEM(index_breccia +1)-delta)/(TI_breccia**2))))372 do iloop =index_breccia + 2,index_bedrock368 TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ & 369 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + & 370 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2)))) 371 do iloop = index_breccia + 2,index_bedrock 373 372 TI_PEM(ig,iloop,islope) = TI_breccia 374 373 enddo 375 374 else ! we keep the high ti values 376 do iloop =index_breccia + 1,index_bedrock375 do iloop = index_breccia + 1,index_bedrock 377 376 TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope) 378 377 enddo … … 380 379 !! transition 381 380 delta = depth_bedrock 382 TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock +1)-layer_PEM(index_bedrock))/&383 (((delta -layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &384 ((layer_PEM(index_bedrock +1)-delta)/(TI_breccia**2))))381 TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ & 382 (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + & 383 ((layer_PEM(index_bedrock + 1) - delta)/(TI_breccia**2)))) 385 384 do iloop = index_bedrock + 2,nsoil_PEM 386 385 TI_PEM(ig,iloop,islope) = TI_bedrock … … 398 397 inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ & 399 398 (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + & 400 ((layer_PEM(index_breccia + 1) -delta)/(TI_breccia**2))))399 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2)))) 401 400 do iloop = index_breccia + 2,index_bedrock 402 401 inertiedat_PEM(ig,iloop) = TI_breccia … … 448 447 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) 449 448 enddo 450 write(*,*) 'PEMETAT0: I ce tabledone'449 write(*,*) 'PEMETAT0: ICE TABLE done' 451 450 else if (icetable_dynamic) then 452 451 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM) … … 454 453 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) 455 454 enddo 456 write(*,*) 'PEMETAT0: I ce tabledone'455 write(*,*) 'PEMETAT0: ICE TABLE done' 457 456 endif 458 457 … … 468 467 endif 469 468 470 write(*,*) 'PEMETAT0: CO2 adsorption done'469 write(*,*) 'PEMETAT0: CO2 & H2O adsorption done' 471 470 endif !soil_pem 472 471 endif ! of if (startphy_file) -
trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90
r3599 r3602 27 27 ! Purpose: Read initial confitions file from the PCM 28 28 ! 29 ! Authors: RV & LL29 ! Authors: JBC 30 30 !======================================================================= 31 31 … … 52 52 real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density_timeseries ! Water density timeseries in the soil layer [kg/m^3] 53 53 ! Local variables 54 integer :: i, j, l, t, islope ! Loop variables 55 real :: A, B, mmean ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer 56 character(2) :: num ! Hor reading sloped variables 57 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen) :: h2o_ice_s_dyn ! H2O ice per slope of the concatenated file [kg/m^2] 58 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen) :: watercap 59 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen) :: perennial_co2ice 60 real, dimension(iim_input + 1,jjm_input + 1,timelen) :: vmr_co2_dyn ! CO2 volume mixing ratio in the first layer [mol/m^3] 61 real, dimension(iim_input + 1,jjm_input + 1,timelen) :: ps_dyn ! Surface pressure [Pa] 62 real, dimension(iim_input + 1,jjm_input + 1) :: ps_avg_dyn ! Averaged surface pressure of the concatenated file [Pa] 63 real, dimension(iim_input + 1,jjm_input + 1,nslope) :: min_co2_ice_dyn 64 real, dimension(iim_input + 1,jjm_input + 1,nslope) :: min_h2o_ice_dyn 65 real, dimension(iim_input + 1,jjm_input + 1,nslope) :: tsurf_avg_dyn ! Averaged surface temperature of the concatenated file [K] 66 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope) :: tsoil_avg_dyn ! Averaged soil temperature of the concatenated file [K] 67 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen) :: tsurf_dyn ! Surface temperature of the concatenated file, time series [K] 68 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_dyn ! Soil temperature of the concatenated file, time series [K] 69 real, dimension(iim_input + 1,jjm_input + 1,timelen) :: q_co2_dyn ! CO2 mass mixing ratio in the first layer [kg/m^3] 70 real, dimension(iim_input + 1,jjm_input + 1,timelen) :: q_h2o_dyn ! H2O mass mixing ratio in the first layer [kg/m^3] 71 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen) :: co2_ice_slope_dyn ! co2 ice amount per slope of the year [kg/m^2] 72 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen) :: watersurf_density_dyn ! Water density at the surface, time series [kg/m^3] 73 real, dimension(iim_input + 1,jjm_input + 1,nslope) :: watersurf_density_dyn_avg ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3] 74 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn ! Water density in the soil layer, time series [kg/m^3] 54 integer :: i, j, l, islope ! Loop variables 55 real :: A, B ! Intermediate variables to compute the mean molar mass of the layer 56 character(2) :: num ! For reading sloped variables 57 58 real, dimension(:,:), allocatable :: var2_read ! Variables for reading (2 dimensions) 59 real, dimension(:,:,:), allocatable :: var3_read_1, var3_read_2, var3_read_3 ! Variables for reading (3 dimensions) 60 real, dimension(:,:,:,:), allocatable :: var4_read ! Variables for reading (4 dimensions) 75 61 !----------------------------------------------------------------------- 76 62 77 !------------------ Year 1------------------63 !------------------------------- Year 1 -------------------------------- 78 64 ! Open the NetCDF file 79 65 write(*,*) "Opening "//filename1//"..." … … 81 67 82 68 ! Download the data from the file 69 allocate(var2_read(iim_input + 1,jjm_input + 1),var3_read_1(iim_input + 1,jjm_input + 1,timelen),var3_read_2(iim_input + 1,jjm_input + 1,timelen)) 70 83 71 if (nslope == 1) then ! There is no slope 84 call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:)) 72 ! CO2 ice 73 !-------- 74 call get_var3("co2ice",var3_read_1) 75 where (var3_read_1 < 0.) var3_read_1 = 0. 85 76 write(*,*) "Data for co2_ice downloaded." 86 87 call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:)) 77 #ifndef CPP_STD 78 call get_var3("perennial_co2ice",var3_read_2) 79 write(*,*) "Data for perennial_co2ice downloaded." 80 #endif 81 82 ! Compute the minimum over the year for each point 83 var2_read = minval(var3_read_1 + var3_read_2,3) 84 #ifndef CPP_1D 85 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,1,1)) 86 #else 87 min_co2_ice(1,1,1) = var2_read(1,1) 88 #endif 89 write(*,*) "Min of co2_ice computed." 90 91 ! H2O ice 92 !-------- 93 call get_var3("h2o_ice_s",var3_read_1) 94 where (var3_read_1 < 0.) var3_read_1 = 0. 88 95 write(*,*) "Data for h2o_ice_s downloaded." 89 90 call get_var3("tsurf",tsurf_dyn(:,:,1,:)) 96 #ifndef CPP_STD 97 call get_var3("watercap",var3_read_2) 98 write(*,*) "Data for watercap downloaded." 99 #endif 100 101 ! Compute the minimum over the year for each point 102 var2_read = minval(var3_read_1 + var3_read_2,3) 103 #ifndef CPP_1D 104 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,1,1)) 105 #else 106 min_h2o_ice(1,1,1) = var2_read(1,1) 107 #endif 108 write(*,*) "Min of h2o_ice computed." 109 110 ! Tsurf 111 !------ 112 call get_var3("tsurf",var3_read_1) 91 113 write(*,*) "Data for tsurf downloaded." 92 114 93 #ifndef CPP_STD 94 call get_var3("watercap",watercap(:,:,1,:)) 95 write(*,*) "Data for watercap downloaded." 96 97 call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:)) 98 write(*,*) "Data for perennial_co2ice downloaded." 99 #endif 115 ! Compute average over the year for each point 116 var2_read = sum(var3_read_1,3)/timelen 117 #ifndef CPP_1D 118 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,1)) 119 #else 120 tsurf_avg_yr1(1,1) = var2_read(1,1) 121 #endif 122 write(*,*) "Average of tsurf computed." 100 123 else ! We use slopes 101 124 do islope = 1,nslope 102 125 write(num,'(i2.2)') islope 103 call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:)) 126 127 ! CO2 ice 128 !-------- 129 call get_var3("co2ice_slope"//num,var3_read_1) 130 where (var3_read_1 < 0.) var3_read_1 = 0. 104 131 write(*,*) "Data for co2_ice_slope"//num//" downloaded." 105 call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:)) 132 #ifndef CPP_STD 133 call get_var3("perennial_co2ice_slope"//num,var3_read_2) 134 write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded." 135 #endif 136 137 ! Compute the minimum over the year for each point 138 var2_read = minval(var3_read_1 + var3_read_2,3) 139 #ifndef CPP_1D 140 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,1)) 141 #else 142 min_co2_ice(1,islope,1) = var2_read(1,1) 143 #endif 144 write(*,*) "Min of co2_ice computed for slope "//num//"." 145 146 ! H2O ice 147 !-------- 148 call get_var3("h2o_ice_s_slope"//num,var3_read_1) 149 where (var3_read_1 < 0.) var3_read_1 = 0. 106 150 write(*,*) "Data for h2o_ice_s_slope"//num//" downloaded." 107 call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:)) 151 #ifndef CPP_STD 152 call get_var3("watercap_slope"//num,var3_read_2) 153 write(*,*) "Data for watercap_slope"//num//" downloaded." 154 #endif 155 156 ! Compute the minimum over the year for each point 157 var2_read = minval(var3_read_1 + var3_read_2,3) 158 #ifndef CPP_1D 159 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,1)) 160 #else 161 min_h2o_ice(1,islope,1) = var2_read(1,1) 162 #endif 163 write(*,*) "Min of h2o_ice computed for slope "//num//"." 164 165 ! Tsurf 166 !------ 167 call get_var3("tsurf_slope"//num,var3_read_1) 108 168 write(*,*) "Data for tsurf_slope"//num//" downloaded." 109 #ifndef CPP_STD 110 call get_var3("watercap_slope"//num,watercap(:,:,islope,:)) 111 write(*,*) "Data for watercap_slope"//num//" downloaded." 112 call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:)) 113 write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded." 114 #endif 169 170 ! Compute average over the year for each point 171 var2_read = sum(var3_read_1,3)/timelen 172 #ifndef CPP_1D 173 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,islope)) 174 #else 175 tsurf_avg_yr1(1,islope) = var2_read(1,1) 176 #endif 177 write(*,*) "Average of tsurf computed for slope "//num//"." 115 178 enddo 116 179 endif 117 180 118 181 ! Close the NetCDF file 119 write(*,*) "Closing "//filename1//"..."120 182 call error_msg(nf90_close(fID),"close",filename1) 121 122 ! Compute the minimum over the year for each point 123 write(*,*) "Computing the min of h2o_ice..." 124 where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0. 125 min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4) 126 write(*,*) "Computing the min of co2_ice..." 127 where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0. 128 min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4) 129 130 ! Compute averages over the year for each point 131 write(*,*) "Computing the average of tsurf..." 132 tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen 133 134 #ifndef CPP_1D 135 call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg) 136 do islope = 1,nslope 137 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope,1)) 138 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,1)) 139 enddo 140 #else 141 tsurf_avg_yr1(1,:) = tsurf_avg_dyn(1,1,:) 142 min_co2_ice(1,:,1) = min_co2_ice_dyn(1,1,:) 143 min_h2o_ice(1,:,1) = min_h2o_ice_dyn(1,1,:) 144 #endif 145 write(*,*) "Downloading data Y1 done!" 146 147 !------------------ Year 2 ------------------ 148 183 write(*,*) '"'//filename1//'" processed!' 184 185 !------------------------------- Year 2 -------------------------------- 149 186 ! Open the NetCDF file 150 187 write(*,*) "Opening "//filename2//"..." … … 152 189 153 190 ! Download the data from the file 154 call get_var3("ps",ps_dyn) 191 allocate(var3_read_3(iim_input + 1,jjm_input + 1,nsoilmx),var4_read(iim_input + 1,jjm_input + 1,nsoilmx,timelen)) 192 193 ! Surface pressure 194 !----------------- 195 call get_var3("ps",var3_read_1) 155 196 write(*,*) "Data for surface pressure downloaded." 156 197 157 call get_var3("co2_layer1",q_co2_dyn) 158 write(*,*) "Data for vmr co2 downloaded." 159 160 call get_var3("h2o_layer1",q_h2o_dyn) 161 write(*,*) "Data for vmr h2o downloaded." 198 ! Compute average over the year for each point 199 var2_read = sum(var3_read_1,3)/timelen 200 #ifndef CPP_1D 201 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,ps_timeseries) 202 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,ps_avg) 203 #else 204 ps_timeseries(1,:) = var3_read_1(1,:) 205 ps_avg(1) = var2_read(1,1) 206 #endif 207 write(*,*) "Average of surface pressure computed." 208 209 ! CO2 vmr 210 !-------- 211 call get_var3("co2_layer1",var3_read_1) 212 where (var3_read_1 < 0.) 213 var3_read_1 = 1.e-10 214 else where (var3_read_1 > 1.) 215 var3_read_1 = 1. 216 end where 217 #ifndef CPP_1D 218 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,q_co2) 219 #else 220 q_co2(1,:) = var3_read_1(1,1,:) 221 #endif 222 A = (1./m_co2 - 1./m_noco2) 223 B = 1./m_noco2 224 vmr_co2_PCM = q_co2/(A*q_co2 + B)/m_co2 225 write(*,*) "Data for CO2 vmr downloaded." 226 227 ! H2O vmr 228 !-------- 229 call get_var3("h2o_layer1",var3_read_1) 230 where (var3_read_1 < 0.) 231 var3_read_1 = 1.e-10 232 else where (var3_read_1 > 1.) 233 var3_read_1 = 1. 234 end where 235 #ifndef CPP_1D 236 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,q_h2o) 237 #else 238 q_h2o(1,:) = var3_read_1(1,1,:) 239 #endif 240 write(*,*) "Data for H2O vmr downloaded." 162 241 163 242 if (nslope == 1) then ! There is no slope 164 call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:)) 243 ! CO2 ice 244 !-------- 245 call get_var3("co2ice",var3_read_1) 246 where (var3_read_1 < 0.) var3_read_1 = 0. 165 247 write(*,*) "Data for co2_ice downloaded." 166 167 call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:)) 248 #ifndef CPP_STD 249 call get_var3("perennial_co2ice",var3_read_2) 250 write(*,*) "Data for perennial_co2ice downloaded." 251 #endif 252 253 ! Compute the minimum over the year for each point 254 var2_read = minval(var3_read_1 + var3_read_2,3) 255 #ifndef CPP_1D 256 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,1,1)) 257 #else 258 min_co2_ice(1,1,1) = var2_read(1,1) 259 #endif 260 write(*,*) "Min of co2_ice computed." 261 262 ! H2O ice 263 !-------- 264 call get_var3("h2o_ice_s",var3_read_1) 265 where (var3_read_1 < 0.) var3_read_1 = 0. 168 266 write(*,*) "Data for h2o_ice_s downloaded." 169 170 call get_var3("tsurf",tsurf_dyn(:,:,1,:)) 267 #ifndef CPP_STD 268 call get_var3("watercap",var3_read_2) 269 write(*,*) "Data for watercap downloaded." 270 #endif 271 272 ! Compute the minimum over the year for each point 273 var2_read = minval(var3_read_1 + var3_read_2,3) 274 #ifndef CPP_1D 275 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,1,1)) 276 #else 277 min_h2o_ice(1,1,1) = var2_read(1,1) 278 #endif 279 write(*,*) "Min of h2o_ice computed." 280 281 ! Tsurf 282 !------ 283 call get_var3("tsurf",var3_read_1) 171 284 write(*,*) "Data for tsurf downloaded." 172 285 173 #ifndef CPP_STD 174 call get_var3("watercap",watercap(:,:,1,:)) 175 write(*,*) "Data for watercap downloaded." 176 177 call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:)) 178 write(*,*) "Data for perennial_co2ice downloaded." 286 ! Compute average over the year for each point 287 var2_read = sum(var3_read_1,3)/timelen 288 #ifndef CPP_1D 289 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,1)) 290 #else 291 tsurf_avg_yr1(1,1) = var2_read(1,1) 292 #endif 293 write(*,*) "Average of tsurf computed." 179 294 180 295 if (soil_pem) then 181 call get_var4("soiltemp",tsoil_dyn(:,:,:,1,:)) 296 ! Tsoil 297 !------ 298 call get_var4("soiltemp",var4_read) 182 299 write(*,*) "Data for soiltemp downloaded." 183 300 184 call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:)) 301 ! Compute average over the year for each point 302 var3_read_3 = sum(var4_read,4)/timelen 303 #ifndef CPP_1D 304 do l = 1,nsoilmx 305 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),tsoil_timeseries(:,l,1,:)) 306 enddo 307 call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var3_read_3,tsoil_avg(:,:,1)) 308 #else 309 tsoil_timeseries(1,:,1,:) = var4_read(1,1,:,:) 310 tsoil_avg(1,:,1) = var3_read_3(1,:,1) 311 #endif 312 write(*,*) "Average of tsoil computed." 313 314 ! Soil water density 315 !------------------- 316 call get_var4("waterdensity_soil",var4_read) 317 #ifndef CPP_1D 318 do l = 1,nsoilmx 319 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),watersoil_density_timeseries(:,l,1,:)) 320 enddo 321 #else 322 watersoil_density_timeseries(1,:,1,:) = var4_read(1,1,:,:) 323 #endif 185 324 write(*,*) "Data for waterdensity_soil downloaded." 186 325 187 call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,1,:)) 326 ! Surface water density 327 !---------------------- 328 call get_var3("waterdensity_surface",var3_read_1) 188 329 write(*,*) "Data for waterdensity_surface downloaded." 189 endif !soil_pem 190 #endif 330 331 ! Compute average over the year for each point 332 var2_read = sum(var3_read_1,3)/timelen 333 #ifndef CPP_1D 334 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,watersurf_density_avg(:,1)) 335 #else 336 watersurf_density_avg(1,1) = var2_read(1,1) 337 #endif 338 write(*,*) "Average of waterdensity_surface computed." 339 endif 191 340 else ! We use slopes 192 341 do islope = 1,nslope 193 342 write(num,'(i2.2)') islope 194 call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:)) 343 344 ! CO2 ice 345 !-------- 346 call get_var3("co2ice_slope"//num,var3_read_1) 347 where (var3_read_1 < 0.) var3_read_1 = 0. 195 348 write(*,*) "Data for co2_ice_slope"//num//" downloaded." 196 call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:)) 349 #ifndef CPP_STD 350 call get_var3("perennial_co2ice_slope"//num,var3_read_2) 351 write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded." 352 #endif 353 354 ! Compute the minimum over the year for each point 355 var2_read = minval(var3_read_1 + var3_read_2,3) 356 #ifndef CPP_1D 357 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,2)) 358 #else 359 min_co2_ice(1,islope,2) = var2_read(1,1) 360 #endif 361 write(*,*) "Min of co2_ice computed for slope "//num//"." 362 363 ! H2O ice 364 !-------- 365 call get_var3("h2o_ice_s_slope"//num,var3_read_1) 366 where (var3_read_1 < 0.) var3_read_1 = 0. 197 367 write(*,*) "Data for h2o_ice_s_slope"//num//" downloaded." 198 call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:)) 368 #ifndef CPP_STD 369 call get_var3("watercap_slope"//num,var3_read_2) 370 write(*,*) "Data for watercap_slope"//num//" downloaded." 371 #endif 372 373 ! Compute the minimum over the year for each point 374 var2_read = minval(var3_read_1 + var3_read_2,3) 375 #ifndef CPP_1D 376 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,2)) 377 #else 378 min_h2o_ice(1,islope,2) = var2_read(1,1) 379 #endif 380 write(*,*) "Min of h2o_ice computed for slope "//num//"." 381 382 ! Tsurf 383 !------ 384 call get_var3("tsurf_slope"//num,var3_read_1) 199 385 write(*,*) "Data for tsurf_slope"//num//" downloaded." 200 #ifndef CPP_STD 201 call get_var3("watercap_slope"//num,watercap(:,:,islope,:)) 202 write(*,*) "Data for watercap_slope"//num//" downloaded." 203 call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:)) 204 write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded." 386 387 ! Compute average over the year for each point 388 var2_read = sum(var3_read_1,3)/timelen 389 #ifndef CPP_1D 390 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg(:,islope)) 391 #else 392 tsurf_avg(1,islope) = var2_read(1,1) 393 #endif 394 write(*,*) "Average of tsurf computed for slope "//num//"." 395 205 396 if (soil_pem) then 206 call get_var4("soiltemp_slope"//num,tsoil_dyn(:,:,:,islope,:)) 397 ! Tsoil 398 !------ 399 call get_var4("soiltemp_slope"//num,var4_read) 207 400 write(*,*) "Data for soiltemp_slope"//num//" downloaded." 208 call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:)) 401 402 ! Compute average over the year for each point 403 var3_read_3 = sum(var4_read,4)/timelen 404 #ifndef CPP_1D 405 do l = 1,nsoilmx 406 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),tsoil_timeseries(:,l,islope,:)) 407 enddo 408 call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var3_read_3,tsoil_avg(:,:,islope)) 409 #else 410 tsoil_timeseries(1,:,islope,:) = var4_read(1,1,:,:) 411 tsoil_avg(1,:,islope) = var3_read_3(1,:,1) 412 #endif 413 write(*,*) "Average of tsoil computed for slope "//num//"." 414 415 ! Soil water density 416 !------------------- 417 call get_var4("waterdensity_soil_slope"//num,var4_read) 418 #ifndef CPP_1D 419 do l = 1,nsoilmx 420 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),watersoil_density_timeseries(:,l,islope,:)) 421 enddo 422 #else 423 watersoil_density_timeseries(1,:,islope,:) = var4_read(1,1,:,:) 424 #endif 209 425 write(*,*) "Data for waterdensity_soil_slope"//num//" downloaded." 210 call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:)) 211 write(*,*) "Data for waterdensity_surface_slope"//num//" downloaded." 426 427 ! Surface water density 428 !---------------------- 429 call get_var3("waterdensity_surface"//num,var3_read_1) 430 write(*,*) "Data for waterdensity_surface"//num//" downloaded." 431 432 ! Compute average over the year for each point 433 var2_read = sum(var3_read_1,3)/timelen 434 #ifndef CPP_1D 435 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,watersurf_density_avg(:,islope)) 436 #else 437 watersurf_density_avg(1,islope) = var2_read(1,1) 438 #endif 439 write(*,*) "Average of waterdensity_surface computed for slope "//num//"." 212 440 endif 213 #endif214 441 enddo 215 442 endif 443 deallocate(var2_read,var3_read_1,var3_read_2,var3_read_3,var4_read) 216 444 217 445 ! Close the NetCDF file 218 write(*,*) "Closing "//filename2//"..."219 446 call error_msg(nf90_close(fID),"close",filename2) 220 221 ! Compute the minimum over the year for each point 222 write(*,*) "Computing the min of h2o_ice..." 223 where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0. 224 min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4) 225 write(*,*) "Computing the min of co2_ice..." 226 where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0. 227 min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4) 228 229 ! Compute averages over the year for each point 230 write(*,*) "Computing the average of tsurf..." 231 tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen 232 233 write(*,*) "Computing the average of Ps..." 234 ps_avg_dyn = sum(ps_dyn,3)/timelen 235 236 if (soil_pem) then 237 write(*,*) "Computing the average of tsoil..." 238 tsoil_avg_dyn = sum(tsoil_dyn,5)/timelen 239 write(*,*) "Computing the average of waterdensity_surface..." 240 watersurf_density_dyn_avg = sum(watersurf_density_dyn,4)/timelen 241 endif 242 243 ! By definition, we get rid of the negative values 244 A = (1./m_co2 - 1./m_noco2) 245 B = 1./m_noco2 246 do i = 1,iim_input + 1 247 do j = 1,jjm_input + 1 248 do t = 1, timelen 249 if (q_co2_dyn(i,j,t) < 0) then 250 q_co2_dyn(i,j,t) = 1.e-10 251 else if (q_co2_dyn(i,j,t) > 1) then 252 q_co2_dyn(i,j,t) = 1. 253 endif 254 if (q_h2o_dyn(i,j,t) < 0) then 255 q_h2o_dyn(i,j,t) = 1.e-10 256 else if (q_h2o_dyn(i,j,t) > 1) then 257 q_h2o_dyn(i,j,t) = 1. 258 endif 259 mmean = 1/(A*q_co2_dyn(i,j,t) + B) 260 vmr_co2_dyn(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2 261 enddo 262 enddo 263 enddo 264 265 #ifndef CPP_1D 266 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_dyn,vmr_co2_PCM) 267 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_dyn,ps_timeseries) 268 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_avg_dyn,ps_avg) 269 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2) 270 call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o) 271 do islope = 1,nslope 272 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope,2)) 273 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,2)) 274 if (soil_pem) then 275 do l = 1,nsoilmx 276 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_avg_dyn(:,:,l,islope),tsoil_avg(:,l,islope)) 277 do t = 1,timelen 278 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_dyn(:,:,l,islope,t),tsoil_timeseries(:,l,islope,t)) 279 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density_timeseries(:,l,islope,t)) 280 enddo 281 enddo 282 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_avg(:,:,islope),watersurf_density_avg(:,islope)) 283 endif ! soil_pem 284 enddo 285 call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg) 286 #else 287 vmr_co2_PCM(1,:) = vmr_co2_dyn(1,1,:) 288 ps_timeseries(1,:) = ps_dyn(1,1,:) 289 ps_avg(1) = ps_avg_dyn(1,1) 290 q_co2(1,:) = q_co2_dyn(1,1,:) 291 q_h2o(1,:) = q_h2o_dyn(1,1,:) 292 min_co2_ice(1,:,2) = min_co2_ice_dyn(1,1,:) 293 min_h2o_ice(1,:,2) = min_h2o_ice_dyn(1,1,:) 294 if (soil_pem) then 295 tsoil_timeseries(1,:,:,:) = tsoil_dyn(1,1,:,:,:) 296 tsoil_avg(1,:,:) = tsoil_avg_dyn(1,1,:,:) 297 watersoil_density_timeseries(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:) 298 watersurf_density_avg(1,:) = watersurf_density_dyn_avg(1,1,:) 299 endif ! soil_pem 300 tsurf_avg(1,:) = tsurf_avg_dyn(1,1,:) 301 #endif 302 write(*,*) "Downloading data Y2 done!" 447 write(*,*) '"'//filename2//'" processed!' 303 448 304 449 END SUBROUTINE read_data_PCM … … 392 537 393 538 END MODULE read_data_PCM_mod 539
Note: See TracChangeset
for help on using the changeset viewer.
