Changeset 3426 for trunk

Timestamp:

Sep 6, 2024, 1:32:48 PM (4 months ago)

Author:

jbclement

Message:

PEM:
Modification in the way of updating the soil temperatures to improve the computation time + Small correction in the calculation of 'timestep'.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• pem.F90 (modified) (9 diffs)

trunk/LMDZ.COMMON/libf/evolution/changelog.txt

@@ -416 +416 @@
 == 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'.

trunk/LMDZ.COMMON/libf/evolution/pem.F90

@@ -204 +204 @@
 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 +341 @@
 year_day = 669
 daysec = 88775.
-timestep = year_day*daysec/year_step
+timestep = year_day*daysec*year_step
 
 !----------------------------- INITIALIZATION --------------------------
@@ -540 +538 @@
 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 +555 @@
 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 +562 @@
 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 +578 @@
 
 if (soil_pem) then
-    allocate(tsoil_anom(ngrid,nsoilmx,nslope))
-    tsoil_anom = tsoil - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - <Tsoil> 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 - <Tsoil> 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 +591 @@
     watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
 endif !soil_pem
-deallocate(tsoil_avg,tsoil_PCM_timeseries)
+deallocate(tsoil_avg)
 
 !------------------------
@@ -895 +894 @@
 
 ! 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 +1073 @@
 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)