Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3608) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3609) @@ -569,2 +569,6 @@ == 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. Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3608) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3609) @@ -84,7 +84,7 @@ 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,6 +406,6 @@ 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,5 +857,4 @@ ! 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,7 +865,7 @@ 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) Index: trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 (revision 3608) +++ trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 (revision 3609) @@ -52,8 +52,7 @@ 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,13 +68,24 @@ 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,18 +93,18 @@ 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,80 +112,24 @@ 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,5 +156,5 @@ 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,13 +194,18 @@ 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,18 +213,18 @@ 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,29 +232,29 @@ 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,143 +262,41 @@ #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