Changeset 3609

Timestamp:

Feb 5, 2025, 2:55:26 PM (6 days ago)

Author:

jbclement

Message:

PEM:

• Simplification of "read_data_PCM_mod.F90" to treat the cases with/without slope.
• Few bug corrections related to 1D.

JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

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

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

@@ -569 +569 @@
 == 03/02/2025 == JBC
 Rescaling the pressure deviation to avoid disproportionate oscillations in case of huge average pressure drop.
+
+== 05/02/2025 == JBC
+- Simplification of "read_data_PCM_mod.F90" to treat the cases with/without slope.
+- Few bug corrections related to 1D. 

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

@@ -84 +84 @@
     use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
     use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
-    use comcstfi_h,         only: mugaz
+    use comcstfi_h,         only: mugaz, r
     use surfini_mod,        only: surfini
-    use comconst_mod,       only: pi, rad, g, r, cpp, kappa
+    use comconst_mod,       only: pi, rad, g, cpp, kappa
 #else
     use tracer_h,           only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
@@ -406 +406 @@
 
     write(*,*) '> Reading "start1D.txt" and "startfi.nc"'
-    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,      &
-                         time_0,pdyn,ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
+    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,         &
+                         time_0,pdyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
                          play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
     nsplit_phys = 1
@@ -857 +857 @@
 ! II_d.3 Update soil temperature
         write(*,*)"> Updating soil temperature profile"
-        allocate(tsoil_avg_old(ngrid,nsoilmx_PEM))
         do islope = 1,nslope
             tsoil_avg_old = tsoil_PEM(:,:,islope)
@@ -866 +865 @@
                 do ig = 1,ngrid
                     do isoil = 1,nsoilmx_PEM
-            ! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
-            tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
-            ! Update of watersoil density
+                        ! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
+                        tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
+                        ! Update of watersoil density
                         watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/tsoil_PEM_timeseries(ig,isoil,islope,t) + alpha_clap_h2o)/tsoil_PEM_timeseries(ig,isoil,islope,t)*mmol(igcm_h2o_vap)/(mugaz*r)
                         if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)

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

@@ -52 +52 @@
 real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density_timeseries ! Water density timeseries in the soil layer [kg/m^3]
 ! Local variables
-integer                                                             :: i, j, l, islope ! Loop variables
-real                                                                :: A, B            ! Intermediate variables to compute the mean molar mass of the layer
-character(2)                                                        :: num             ! For reading sloped variables
-
+integer                               :: i, j, l, islope                       ! Loop variables
+real                                  :: A, B                                  ! Intermediate variables to compute the mean molar mass of the layer
+character(:), allocatable             :: num                                   ! For reading sloped variables
 real, dimension(:,:),     allocatable :: var2_read                             ! Variables for reading (2 dimensions)
 real, dimension(:,:,:),   allocatable :: var3_read_1, var3_read_2, var3_read_3 ! Variables for reading (3 dimensions)
@@ -69 +68 @@
 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))
 
-if (nslope == 1) then ! There is no slope
+if (nslope == 1) then ! No slope
+    allocate(character(0) :: num)
+else ! We use slopes
+    allocate(character(8) :: num)
+endif
+
+do islope = 1,nslope
+    if (nslope /= 1) then
+        write(num,'(i2.2)') islope
+        num = '_slope'//num
+    endif
+
     ! CO2 ice
     !--------
-    call get_var3("co2ice",var3_read_1)
+    call get_var3("co2ice"//num,var3_read_1)
     where (var3_read_1 < 0.) var3_read_1 = 0.
-    write(*,*) "Data for co2_ice downloaded."
+    write(*,*) "Data for co2_ice"//num//" downloaded."
 #ifndef CPP_STD
-    call get_var3("perennial_co2ice",var3_read_2)
-    write(*,*) "Data for perennial_co2ice downloaded."
+    call get_var3("perennial_co2ice"//num,var3_read_2)
+    write(*,*) "Data for perennial_co2ice"//num//" downloaded."
 #endif
 
@@ -83 +93 @@
     var2_read = minval(var3_read_1 + var3_read_2,3)
 #ifndef CPP_1D
-    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,1,1))
-#else
-    min_co2_ice(1,1,1) = var2_read(1,1)
-#endif
-    write(*,*) "Min of co2_ice computed."
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,1))
+#else
+    min_co2_ice(1,islope,1) = var2_read(1,1)
+#endif
+    write(*,*) "Min of co2_ice"//num//" computed."
 
     ! H2O ice
     !--------
-    call get_var3("h2o_ice_s",var3_read_1)
+    call get_var3("h2o_ice_s"//num,var3_read_1)
     where (var3_read_1 < 0.) var3_read_1 = 0.
-    write(*,*) "Data for h2o_ice_s downloaded."
+    write(*,*) "Data for h2o_ice_s"//num//" downloaded."
 #ifndef CPP_STD
-    call get_var3("watercap",var3_read_2)
-    write(*,*) "Data for watercap downloaded."
+    call get_var3("watercap"//num,var3_read_2)
+    write(*,*) "Data for watercap"//num//" downloaded."
 #endif
 
@@ -102 +112 @@
     var2_read = minval(var3_read_1 + var3_read_2,3)
 #ifndef CPP_1D
-    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,1,1))
-#else
-    min_h2o_ice(1,1,1) = var2_read(1,1)
-#endif
-    write(*,*) "Min of h2o_ice computed."
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,1))
+#else
+    min_h2o_ice(1,islope,1) = var2_read(1,1)
+#endif
+    write(*,*) "Min of h2o_ice"//num//" computed."
 
     ! Tsurf
     !------
-    call get_var3("tsurf",var3_read_1)
-    write(*,*) "Data for tsurf downloaded."
+    call get_var3("tsurf"//num,var3_read_1)
+    write(*,*) "Data for tsurf"//num//" downloaded."
 
     ! Compute average over the year for each point
     var2_read = sum(var3_read_1,3)/timelen
 #ifndef CPP_1D
-    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,1))
-#else
-    tsurf_avg_yr1(1,1) = var2_read(1,1)
-#endif
-    write(*,*) "Average of tsurf computed."
-else ! We use slopes
-    do islope = 1,nslope
-        write(num,'(i2.2)') islope
-
-        ! CO2 ice
-        !--------
-        call get_var3("co2ice_slope"//num,var3_read_1)
-        where (var3_read_1 < 0.) var3_read_1 = 0.
-        write(*,*) "Data for co2_ice_slope"//num//" downloaded."
-#ifndef CPP_STD
-        call get_var3("perennial_co2ice_slope"//num,var3_read_2)
-        write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded."
-#endif
-
-        ! Compute the minimum over the year for each point
-        var2_read = minval(var3_read_1 + var3_read_2,3)
-#ifndef CPP_1D
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,1))
-#else
-        min_co2_ice(1,islope,1) = var2_read(1,1)
-#endif
-        write(*,*) "Min of co2_ice computed for slope "//num//"."
-
-        ! H2O ice
-        !--------
-        call get_var3("h2o_ice_s_slope"//num,var3_read_1)
-        where (var3_read_1 < 0.) var3_read_1 = 0.
-        write(*,*) "Data for h2o_ice_s_slope"//num//" downloaded."
-#ifndef CPP_STD
-        call get_var3("watercap_slope"//num,var3_read_2)
-        write(*,*) "Data for watercap_slope"//num//" downloaded."
-#endif
-
-        ! Compute the minimum over the year for each point
-        var2_read = minval(var3_read_1 + var3_read_2,3)
-#ifndef CPP_1D
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,1))
-#else
-        min_h2o_ice(1,islope,1) = var2_read(1,1)
-#endif
-        write(*,*) "Min of h2o_ice computed for slope "//num//"."
-
-        ! Tsurf
-        !------
-        call get_var3("tsurf_slope"//num,var3_read_1)
-        write(*,*) "Data for tsurf_slope"//num//" downloaded."
-
-        ! Compute average over the year for each point
-        var2_read = sum(var3_read_1,3)/timelen
-#ifndef CPP_1D
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,islope))
-#else
-        tsurf_avg_yr1(1,islope) = var2_read(1,1)
-#endif
-        write(*,*) "Average of tsurf computed for slope "//num//"."
-    enddo
-endif
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,islope))
+#else
+    tsurf_avg_yr1(1,islope) = var2_read(1,1)
+#endif
+    write(*,*) "Average of tsurf"//num//" computed."
+enddo
 
 ! Close the NetCDF file
@@ -202 +156 @@
     call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,ps_avg)
 #else
-    ps_timeseries(1,:) = var3_read_1(1,:)
+    ps_timeseries(1,:) = var3_read_1(1,1,:)
     ps_avg(1) = var2_read(1,1)
 #endif
@@ -240 +194 @@
 write(*,*) "Data for H2O vmr downloaded."
 
-if (nslope == 1) then ! There is no slope
+do islope = 1,nslope
+    if (nslope /= 1) then
+        write(num,'(i2.2)') islope
+        num = '_slope'//num
+    endif
+
     ! CO2 ice
     !--------
-    call get_var3("co2ice",var3_read_1)
+    call get_var3("co2ice"//num,var3_read_1)
     where (var3_read_1 < 0.) var3_read_1 = 0.
-    write(*,*) "Data for co2_ice downloaded."
+    write(*,*) "Data for co2_ice"//num//" downloaded."
 #ifndef CPP_STD
-    call get_var3("perennial_co2ice",var3_read_2)
-    write(*,*) "Data for perennial_co2ice downloaded."
+    call get_var3("perennial_co2ice"//num,var3_read_2)
+    write(*,*) "Data for perennial_co2ice"//num//" downloaded."
 #endif
 
@@ -254 +213 @@
     var2_read = minval(var3_read_1 + var3_read_2,3)
 #ifndef CPP_1D
-    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,1,1))
-#else
-    min_co2_ice(1,1,1) = var2_read(1,1)
-#endif
-    write(*,*) "Min of co2_ice computed."
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,2))
+#else
+    min_co2_ice(1,islope,2) = var2_read(1,1)
+#endif
+    write(*,*) "Min of co2_ice"//num//" computed."
 
     ! H2O ice
     !--------
-    call get_var3("h2o_ice_s",var3_read_1)
+    call get_var3("h2o_ice_s"//num,var3_read_1)
     where (var3_read_1 < 0.) var3_read_1 = 0.
-    write(*,*) "Data for h2o_ice_s downloaded."
+    write(*,*) "Data for h2o_ice_s"//num//" downloaded."
 #ifndef CPP_STD
-    call get_var3("watercap",var3_read_2)
-    write(*,*) "Data for watercap downloaded."
+    call get_var3("watercap"//num,var3_read_2)
+    write(*,*) "Data for watercap"//num//" downloaded."
 #endif
 
@@ -273 +232 @@
     var2_read = minval(var3_read_1 + var3_read_2,3)
 #ifndef CPP_1D
-    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,1,1))
-#else
-    min_h2o_ice(1,1,1) = var2_read(1,1)
-#endif
-    write(*,*) "Min of h2o_ice computed."
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,2))
+#else
+    min_h2o_ice(1,islope,2) = var2_read(1,1)
+#endif
+    write(*,*) "Min of h2o_ice"//num//" computed."
 
     ! Tsurf
     !------
-    call get_var3("tsurf",var3_read_1)
-    write(*,*) "Data for tsurf downloaded."
+    call get_var3("tsurf"//num,var3_read_1)
+       write(*,*) "Data for tsurf"//num//" downloaded."
 
     ! Compute average over the year for each point
     var2_read = sum(var3_read_1,3)/timelen
 #ifndef CPP_1D
-    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,1))
-#else
-    tsurf_avg_yr1(1,1) = var2_read(1,1)
-#endif
-    write(*,*) "Average of tsurf computed."
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg(:,islope))
+#else
+    tsurf_avg(1,islope) = var2_read(1,1)
+#endif
+    write(*,*) "Average of tsurf"//num//" computed."
 
     if (soil_pem) then
         ! Tsoil
         !------
-        call get_var4("soiltemp",var4_read)
-        write(*,*) "Data for soiltemp downloaded."
+        call get_var4("soiltemp"//num,var4_read)
+        write(*,*) "Data for soiltemp"//num//" downloaded."
 
         ! Compute average over the year for each point
@@ -303 +262 @@
 #ifndef CPP_1D
         do l = 1,nsoilmx
-            call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),tsoil_timeseries(:,l,1,:))
+            call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),tsoil_timeseries(:,l,islope,:))
         enddo
-        call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var3_read_3,tsoil_avg(:,:,1))
-#else
-        tsoil_timeseries(1,:,1,:) = var4_read(1,1,:,:)
-        tsoil_avg(1,:,1) = var3_read_3(1,:,1)
-#endif
-        write(*,*) "Average of tsoil computed."
+        call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var3_read_3,tsoil_avg(:,:,islope))
+#else
+        tsoil_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
+        tsoil_avg(1,:,islope) = var3_read_3(1,1,:)
+#endif
+        write(*,*) "Average of tsoil"//num//" computed."
 
         ! Soil water density
         !-------------------
-        call get_var4("waterdensity_soil",var4_read)
+        call get_var4("waterdensity_soil"//num,var4_read)
 #ifndef CPP_1D
         do l = 1,nsoilmx
-            call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),watersoil_density_timeseries(:,l,1,:))
+            call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),watersoil_density_timeseries(:,l,islope,:))
         enddo
 #else
-        watersoil_density_timeseries(1,:,1,:) = var4_read(1,1,:,:)
-#endif
-        write(*,*) "Data for waterdensity_soil downloaded."
+        watersoil_density_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
+#endif
+        write(*,*) "Data for waterdensity_soil"//num//" downloaded."
 
         ! Surface water density
         !----------------------
-        call get_var3("waterdensity_surface",var3_read_1)
-        write(*,*) "Data for waterdensity_surface downloaded."
+        call get_var3("waterdensity_surface"//num,var3_read_1)
+        write(*,*) "Data for waterdensity_surface"//num//" downloaded."
 
         ! Compute average over the year for each point
         var2_read = sum(var3_read_1,3)/timelen
 #ifndef CPP_1D
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,watersurf_density_avg(:,1))
-#else
-        watersurf_density_avg(1,1) = var2_read(1,1)
-#endif
-        write(*,*) "Average of waterdensity_surface computed."
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,watersurf_density_avg(:,islope))
+#else
+        watersurf_density_avg(1,islope) = var2_read(1,1)
+#endif
+        write(*,*) "Average of waterdensity_surface"//num//" computed."
     endif
-else ! We use slopes
-    do islope = 1,nslope
-        write(num,'(i2.2)') islope
-
-        ! CO2 ice
-        !--------
-        call get_var3("co2ice_slope"//num,var3_read_1)
-        where (var3_read_1 < 0.) var3_read_1 = 0.
-        write(*,*) "Data for co2_ice_slope"//num//" downloaded."
-#ifndef CPP_STD
-        call get_var3("perennial_co2ice_slope"//num,var3_read_2)
-        write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded."
-#endif
-
-        ! Compute the minimum over the year for each point
-        var2_read = minval(var3_read_1 + var3_read_2,3)
-#ifndef CPP_1D
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,2))
-#else
-        min_co2_ice(1,islope,2) = var2_read(1,1)
-#endif
-        write(*,*) "Min of co2_ice computed for slope "//num//"."
-
-        ! H2O ice
-        !--------
-        call get_var3("h2o_ice_s_slope"//num,var3_read_1)
-        where (var3_read_1 < 0.) var3_read_1 = 0.
-        write(*,*) "Data for h2o_ice_s_slope"//num//" downloaded."
-#ifndef CPP_STD
-        call get_var3("watercap_slope"//num,var3_read_2)
-        write(*,*) "Data for watercap_slope"//num//" downloaded."
-#endif
-
-        ! Compute the minimum over the year for each point
-        var2_read = minval(var3_read_1 + var3_read_2,3)
-#ifndef CPP_1D
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,2))
-#else
-        min_h2o_ice(1,islope,2) = var2_read(1,1)
-#endif
-        write(*,*) "Min of h2o_ice computed for slope "//num//"."
-
-        ! Tsurf
-        !------
-        call get_var3("tsurf_slope"//num,var3_read_1)
-        write(*,*) "Data for tsurf_slope"//num//" downloaded."
-
-        ! Compute average over the year for each point
-        var2_read = sum(var3_read_1,3)/timelen
-#ifndef CPP_1D
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg(:,islope))
-#else
-        tsurf_avg(1,islope) = var2_read(1,1)
-#endif
-        write(*,*) "Average of tsurf computed for slope "//num//"."
-
-        if (soil_pem) then
-            ! Tsoil
-            !------
-            call get_var4("soiltemp_slope"//num,var4_read)
-            write(*,*) "Data for soiltemp_slope"//num//" downloaded."
-
-            ! Compute average over the year for each point
-            var3_read_3 = sum(var4_read,4)/timelen
-#ifndef CPP_1D
-            do l = 1,nsoilmx
-                call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),tsoil_timeseries(:,l,islope,:))
-            enddo
-            call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var3_read_3,tsoil_avg(:,:,islope))
-#else
-            tsoil_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
-            tsoil_avg(1,:,islope) = var3_read_3(1,:,1)
-#endif
-            write(*,*) "Average of tsoil computed for slope "//num//"."
-
-            ! Soil water density
-            !-------------------
-            call get_var4("waterdensity_soil_slope"//num,var4_read)
-#ifndef CPP_1D
-            do l = 1,nsoilmx
-                call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),watersoil_density_timeseries(:,l,islope,:))
-            enddo
-#else
-            watersoil_density_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
-#endif
-            write(*,*) "Data for waterdensity_soil_slope"//num//" downloaded."
-
-            ! Surface water density
-            !----------------------
-            call get_var3("waterdensity_surface"//num,var3_read_1)
-            write(*,*) "Data for waterdensity_surface"//num//" downloaded."
-
-            ! Compute average over the year for each point
-            var2_read = sum(var3_read_1,3)/timelen
-#ifndef CPP_1D
-            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,watersurf_density_avg(:,islope))
-#else
-            watersurf_density_avg(1,islope) = var2_read(1,1)
-#endif
-            write(*,*) "Average of waterdensity_surface computed for slope "//num//"."
-        endif
-    enddo
-endif
-deallocate(var2_read,var3_read_1,var3_read_2,var3_read_3,var4_read)
+enddo
+deallocate(var2_read,var3_read_1,var3_read_2,var3_read_3,var4_read,num)
 
 ! Close the NetCDF file