Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3424) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3426) @@ -416,2 +416,5 @@ == 02/09/2024 == JBC Modification to the layering algorithm to make dust lag prevent further sublimation. + +== 06/09/2024 == JBC +Modification in the way of updating the soil temperatures to improve the computation time + Small correction in the calculation of 'timestep'. Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3424) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3426) @@ -204,10 +204,8 @@ real, dimension(:,:), allocatable :: tsurf_avg ! Physic x SLOPE field: Averaged Surface Temperature [K] real, dimension(:,:,:), allocatable :: tsoil_avg ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K] -real, dimension(:,:,:), allocatable :: tsoil_anom ! Amplitude between instataneous and yearly average soil temperature [K] real, dimension(:,:,:), allocatable :: tsurf_PCM_timeseries ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K] -real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SLOPE x TIMES field: Non averaged Soil Temperature [K] -real, dimension(:,:,:,:), allocatable :: tsoil_PCM_timeseries ! IG x SLOPE x TIMES field: Non averaged Soil Temperature [K] +real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K] +real, dimension(:,:,:,:), allocatable :: tsoil_anom ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K] real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K] -real, dimension(:,:), allocatable :: TI_locslope ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: Intermediate when computing Tsoil [K] real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] @@ -343,5 +341,5 @@ year_day = 669 daysec = 88775. -timestep = year_day*daysec/year_step +timestep = year_day*daysec*year_step !----------------------------- INITIALIZATION -------------------------- @@ -540,5 +538,5 @@ allocate(tsurf_avg(ngrid,nslope)) allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen)) -allocate(tsoil_PCM_timeseries(ngrid,nsoilmx,nslope,timelen)) +allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen)) allocate(q_co2_PEM_phys(ngrid,timelen)) allocate(q_h2o_PEM_phys(ngrid,timelen)) @@ -557,5 +555,5 @@ write(*,*) "Downloading data Y1..." call read_data_PCM("data_PCM_Y1.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,1),min_h2o_ice(:,:,1), & - tsurf_avg_yr1,tsoil_avg,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & + tsurf_avg_yr1,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, & watersurf_density_avg,watersoil_density_timeseries) write(*,*) "Downloading data Y1 done!" @@ -564,5 +562,5 @@ write(*,*) "Downloading data Y2..." call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,2),min_h2o_ice(:,:,2), & - tsurf_avg,tsoil_avg,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & + tsurf_avg,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, & watersurf_density_avg,watersoil_density_timeseries) write(*,*) "Downloading data Y2 done!" @@ -580,9 +578,10 @@ if (soil_pem) then - allocate(tsoil_anom(ngrid,nsoilmx,nslope)) - tsoil_anom = tsoil - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - to recompute properly tsoil in the restart + do t = 1,timelen + tsoil_anom(:,:,:,t) = tsoil_anom(:,:,:,t) - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - to recompute properly tsoil in the restart + enddo call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg - tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_PCM_timeseries + tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_anom watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries do l = nsoilmx + 1,nsoilmx_PEM @@ -592,5 +591,5 @@ watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen endif !soil_pem -deallocate(tsoil_avg,tsoil_PCM_timeseries) +deallocate(tsoil_avg) !------------------------ @@ -895,32 +894,26 @@ ! II_d.2 Update soil temperature - allocate(TI_locslope(ngrid,nsoilmx_PEM)) + write(*,*)"Updating soil temperature" allocate(Tsoil_locslope(ngrid,nsoilmx_PEM)) - allocate(Tsurf_locslope(ngrid)) - write(*,*)"Updating soil temperature" - - ! Soil averaged do islope = 1,nslope - TI_locslope = TI_PEM(:,:,islope) + call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) + call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) + do t = 1,timelen - Tsoil_locslope = tsoil_phys_PEM_timeseries(:,:,islope,t) - Tsurf_locslope = tsurf_PCM_timeseries(:,islope,t) - call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) - call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) - tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope + Tsoil_locslope(:,1:nsoilmx) = tsoil_PEM(:,1:nsoilmx,islope) + tsoil_anom(:,:,islope,t) + Tsoil_locslope(:,nsoilmx + 1:) = tsoil_PEM(:,nsoilmx + 1:,islope) do ig = 1,ngrid do isoil = 1,nsoilmx_PEM 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) - if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1) + if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1) enddo enddo enddo enddo - tsoil_PEM = sum(tsoil_phys_PEM_timeseries,4)/timelen watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen write(*,*) "Update of soil temperature done" - deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope) + deallocate(Tsoil_locslope) ! II_d.3 Update the ice table @@ -1080,10 +1073,10 @@ if (soil_pem) then call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) - tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom - deallocate(tsoil_anom) + tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen) #ifndef CPP_STD flux_geo = fluxgeo #endif endif +deallocate(tsoil_anom) ! III_a.4 Pressure (for start)