Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3602) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3603) @@ -562,2 +562,6 @@ - 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. - Few cleanings. + +== 28/01/2025 == JBC +- Bug correction about the pressure/teta reconstruction for the PCM (mismatch between the physical and dynamical grid). +- Improvement of messages giving info about the PEM workflow in the terminal. Index: trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90 (revision 3602) +++ trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90 (revision 3603) @@ -34,5 +34,5 @@ !======================================================================= ! Evolution of CO2 ice for each physical point -write(*,*) 'Evolution of co2 ice' +write(*,*) '> Evolution of CO2 ice' co2_ice_old = co2_ice @@ -88,4 +88,5 @@ real, dimension(ngrid,nslope) :: new_tend ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done !======================================================================= +write(*,*) '> Evolution of H2O ice' if (ngrid /= 1) then ! Not in 1D ! We compute the amount of condensing and sublimating h2o ice Index: trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90 (revision 3602) +++ trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90 (revision 3603) @@ -53,5 +53,5 @@ real, dimension(ngrid,nslope) :: hmax ! Grid points x Slope field: maximum thickness for co2 glacier before flow -write(*,*) "Flow of CO2 glaciers" +write(*,*) "> Flow of CO2 glaciers" call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond) call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax) @@ -93,5 +93,5 @@ real, dimension(ngrid,nslope) :: hmax ! Grid points x Slope field: maximum thickness for co2 glacier before flow -write(*,*) "Flow of H2O glaciers" +write(*,*) "> Flow of H2O glaciers" call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax) call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,h2oice,flag_h2oflow) Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3602) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3603) @@ -139,5 +139,5 @@ real, dimension(ip1jmp1,llm) :: teta ! Potential temperature real, dimension(:,:,:), allocatable :: q ! champs advectes -real, dimension(ip1jmp1) :: ps_start_dyn ! surface pressure in the start file (dynamic grid) +real, dimension(:), allocatable :: pdyn ! pressure for the dynamic grid real, dimension(:), allocatable :: ps_start ! surface pressure in the start file real, dimension(:), allocatable :: ps_start0 ! surface pressure in the start file at the beginning @@ -286,4 +286,13 @@ integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap logical :: num_str + +write(*,*) ' * . . + . * .. + . . . ' +write(*,*) ' + _______ ________ ____ ____ * + ' +write(*,*) ' + . * |_ __ \|_ __ ||_ \ / _| . *' +write(*,*) ' . . | |__) | | |_ \_| | \/ | * * . ' +write(*,*) ' . | ___/ | _| _ | |\ /| | . . ' +write(*,*) '. * * _| |_ _| |__/ | _| |_\/_| |_ * ' +write(*,*) ' |_____| |________||_____||_____| + . ' +write(*,*) ' . * . * .. + * . + .' ! Elapsed time with system clock @@ -371,4 +380,5 @@ ! I_a Read the "run.def" !------------------------ +write(*,*) '> Reading "run.def" (PEM)' #ifndef CPP_1D dtphys = daysec/48. ! Dummy value (overwritten in phyetat0) @@ -379,5 +389,5 @@ allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid)) #else - allocate(q(1,llm,nqtot)) + allocate(q(1,llm,nqtot),pdyn(1)) allocate(longitude(1),latitude(1),cell_area(1)) @@ -395,8 +405,8 @@ endif - call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, & - time_0,ps_start_dyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, & + 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, & play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau) - ps_start_dyn(2) = ps_start_dyn(1) nsplit_phys = 1 day_step = steps_per_sol @@ -410,9 +420,11 @@ !------------------------ ! I_b.1 Read "start.nc" +write(*,*) '> Reading "start.nc"' allocate(ps_start0(ngrid)) #ifndef CPP_1D - call dynetat0(start_name,vcov,ucov,teta,q,masse,ps_start_dyn,phis,time_0) - - call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps_start_dyn,ps_start0) + allocate(pdyn(ip1jmp1)) + call dynetat0(start_name,vcov,ucov,teta,q,masse,pdyn,phis,time_0) + + call gr_dyn_fi(1,iip1,jjp1,ngridmx,pdyn,ps_start0) call iniconst ! Initialization of dynamical constants (comconst_mod) @@ -431,6 +443,7 @@ call iniphysiq(iim,jjm,llm,(jjm - 1)*iim + 2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys) #else - ps_start0(1) = ps_start_dyn(1) + ps_start0 = pdyn #endif +deallocate(pdyn) ! In the PCM, these values are given to the physic by the dynamic. @@ -450,4 +463,5 @@ ! First we read the initial state (starfi.nc) #ifndef CPP_STD + write(*,*) '> Reading "startfi.nc"' call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, & tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, & @@ -510,5 +524,5 @@ #endif -do nnq = 1,nqtot ! Why not using ini_tracer ? +do nnq = 1,nqtot ! Why not using ini_tracer? if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq if (noms(nnq) == "h2o_vap") then @@ -529,5 +543,5 @@ enddo write(*,*) 'Flat slope for islope = ',iflat -write(*,*) 'corresponding criterium = ',def_slope_mean(iflat) +write(*,*) 'Corresponding criterium = ',def_slope_mean(iflat) !------------------------ @@ -597,5 +611,5 @@ ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface ps_avg_global_new = ps_avg_global_ini -write(*,*) "Total surface of the planet =", total_surface +write(*,*) "Total surface of the planet =", total_surface write(*,*) "Initial global averaged pressure =", ps_avg_global_ini @@ -604,4 +618,5 @@ ! I_h Read the "startpem.nc" !------------------------ +write(*,*) '> Reading "startpem.nc"' allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),stratif(ngrid,nslope)) allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid)) @@ -642,7 +657,6 @@ enddo enddo -write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_sublim_surf_ini -write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf - +write(*,*) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini +write(*,*) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf if (adsorption_pem) then @@ -707,5 +721,5 @@ do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear) ! II.a.1. Compute updated global pressure - write(*,*) "Recomputing the new pressure..." + write(*,*) "> Updating the surface pressure" ps_avg_global_old = ps_avg_global_new do i = 1,ngrid @@ -720,10 +734,10 @@ enddo endif - write(*,*) 'Global average pressure old time step',ps_avg_global_old - write(*,*) 'Global average pressure new time step',ps_avg_global_new + write(*,*) 'Global average pressure old time step:',ps_avg_global_old + write(*,*) 'Global average pressure new time step:',ps_avg_global_new ! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption) + write(*,*) "> Updating the surface pressure timeseries for the new pressure" allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen)) - write(*,*) "Recomputing the pressure levels timeseries adapted to the old pressure..." do l = 1,nlayer + 1 do ig = 1,ngrid @@ -731,9 +745,8 @@ enddo enddo - write(*,*) "Recomputing the surface pressure timeseries adapted to the new pressure..." do ig = 1,ngrid ps_timeseries(ig,:) = ps_timeseries(ig,:)*ps_avg_global_new/ps_avg_global_old enddo - write(*,*) "Recomputing the pressure levels timeseries adapted to the new pressure..." + write(*,*) "> Updating the pressure levels timeseries for the new pressure" allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen)) do l = 1,nlayer + 1 @@ -744,5 +757,5 @@ ! II.a.3. Tracers timeseries - write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..." + write(*,*) "> Updating the tracer VMR timeseries for the new pressure" allocate(vmr_co2_PEM_phys(ngrid,timelen)) l = 1 @@ -805,5 +818,5 @@ !------------------------ ! II_d.1 Update Tsurf - write(*,*) "Updating the new Tsurf" + write(*,*) "> Updating surface temperature" do ig = 1,ngrid do islope = 1,nslope @@ -811,5 +824,5 @@ if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then co2ice_disappeared(ig,islope) = .true. - if (latitude_deg(ig) > 0) then + if (latitude_deg(ig) > 0.) then outer1: do ig_loop = ig,ngrid do islope_loop = islope,iflat,-1 @@ -838,9 +851,10 @@ if (soil_pem) then ! II_d.2 Shifting soil temperature to surface + write(*,*) "> Shifting soil temperature profile to match surface evolution" call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM) deallocate(zshift_surf,zlag) ! II_d.3 Update soil temperature - write(*,*)"Updating soil temperature" + write(*,*)"> Updating soil temperature profile" allocate(tsoil_avg_old(ngrid,nsoilmx_PEM)) do islope = 1,nslope @@ -863,15 +877,14 @@ watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen deallocate(tsoil_avg_old) - write(*,*) "Update of soil temperature done" ! II_d.4 Update the ice table allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope)) if (icetable_equilibrium) then - write(*,*) "Compute ice table (equilibrium method)" + write(*,*) "> Updating ice table (equilibrium method)" icetable_thickness_old = icetable_thickness call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness) call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere else if (icetable_dynamic) then - write(*,*) "Compute ice table (dynamic method)" + write(*,*) "> Updating ice table (dynamic method)" ice_porefilling_old = ice_porefilling allocate(porefill(nsoilmx_PEM)) @@ -916,6 +929,6 @@ enddo enddo - write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbed - write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbed + write(*,*) "Total mass of CO2 in the regolith =", totmassco2_adsorbed + write(*,*) "Total mass of H2O in the regolith =", totmassh2o_adsorbed endif endif !soil_pem @@ -965,5 +978,5 @@ ! II_g Checking the stopping criterion !------------------------ - write(*,*) "Checking the stopping criteria..." + write(*,*) "> Checking the stopping criteria" call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid) call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope) @@ -998,6 +1011,6 @@ else write(*,*) "The PEM can continue!" - write(*,*) "Number of iterations of the PEM: i_myear_leg =", i_myear_leg - write(*,*) "Number of simulated Martian years: i_myear =", i_myear + write(*,*) "Number of iterations of the PEM : i_myear_leg =", i_myear_leg + write(*,*) "Number of simulated Martian years: i_myear =", i_myear endif enddo ! big time iteration loop of the pem @@ -1024,4 +1037,5 @@ !------------------------ ! III_a.1 Ice update for start file +write(*,*) '> Reconstructing perennial ice and frost for the PCM' watercap = 0. perennial_co2ice = co2_ice @@ -1052,4 +1066,5 @@ ! III.a.2. Tsurf update for start file +write(*,*) '> Reconstructing the surface temperature for the PCM' tsurf = tsurf_avg + tsurf_dev deallocate(tsurf_dev) @@ -1057,4 +1072,5 @@ ! III_a.3 Tsoil update for start file if (soil_pem) then + write(*,*) '> Reconstructing the soil temperature profile for the PCM' inertiesoil = TI_PEM(:,:nsoilmx,:) ! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value @@ -1070,4 +1086,5 @@ ! III_a.4 Pressure update for start file +write(*,*) '> Reconstructing the pressure for the PCM' allocate(ps_start(ngrid)) ps_start = ps_avg + ps_dev @@ -1075,4 +1092,5 @@ ! III_a.5 Tracers update for start file +write(*,*) '> Reconstructing the tracer VMR for the PCM' allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1)) do l = 1,nlayer + 1 @@ -1118,4 +1136,5 @@ ! III_a.6 Albedo update for start file +write(*,*) '> Reconstructing the albedo for the PCM' do ig = 1,ngrid if (latitude(ig) < 0.) then @@ -1124,5 +1143,5 @@ icap = 1 ! Northern hemisphere endif - do islope = 1,ngrid + do islope = 1,nslope ! Bare ground albedo(ig,:,islope) = albedodat(ig) @@ -1154,4 +1173,5 @@ ! III_a.7 Orbital parameters update for start file +write(*,*) '> Setting the new orbital parameters' if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg) @@ -1164,28 +1184,29 @@ pday = day_ini ztime_fin = time_phys - #ifndef CPP_1D + write(*,*) '> Writing "restart.nc"' + ! Correction on teta due to surface pressure changes + allocate(pdyn(ip1jmp1)) + call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start0/ps_start,pdyn) + do i = 1,ip1jmp1 + teta(i,:) = teta(i,:)*pdyn(i)**kappa + enddo + ! Correction on atmospheric pressure allocate(p(ip1jmp1,nlayer + 1)) - ! Correction on teta due to surface pressure changes - do l = 1,nlayer - do i = 1,ip1jmp1 - teta(i,l)= teta(i,l)*(ps_start0(i)/ps_start(i))**kappa - enddo - enddo - ! Correction on atmospheric pressure - call pression(ip1jmp1,ap,bp,ps_start,p) + call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start,pdyn) + call pression(ip1jmp1,ap,bp,pdyn,p) ! Correction on the mass of atmosphere call massdair(p,masse) call dynredem0("restart.nc",day_ini,phis) call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps_start) - write(*,*) "restart.nc has been written." - deallocate(ap,bp,p) + deallocate(ap,bp,p,pdyn) #else + write(*,*) '> Writing "restart1D.txt"' call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q) - write(*,*) "restart1D.txt has been written." #endif deallocate(ps_start0,ps_start) ! III_b.2 Write the "restartfi.nc" +write(*,*) '> Writing "restartfi.nc"' #ifndef CPP_STD call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, & @@ -1207,5 +1228,4 @@ tsea_ice,sea_ice) #endif -write(*,*) "restartfi.nc has been written." !------------------------ @@ -1213,16 +1233,16 @@ ! III_c Write the "restartpem.nc" !------------------------ +write(*,*) '> Writing "restartpem.nc"' if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings) call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist) call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, & co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,stratif) -write(*,*) "restartpem.nc has been written." call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear) -write(*,*) "The PEM has run for", i_myear_leg, "Martian years." -write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years." -write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years." -write(*,*) "PEM: so far, so good!" +write(*,*) ">> The PEM has run for", i_myear_leg, "Martian years." +write(*,*) ">> The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years." +write(*,*) ">> The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years." +write(*,*) ">> PEM: so far, so good!" if (layering_algo) then Index: trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 (revision 3602) +++ trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 (revision 3603) @@ -181,5 +181,5 @@ ! Close the NetCDF file call error_msg(nf90_close(fID),"close",filename1) -write(*,*) '"'//filename1//'" processed!' +write(*,*) '> "'//filename1//'" processed!' !------------------------------- Year 2 -------------------------------- @@ -445,5 +445,5 @@ ! Close the NetCDF file call error_msg(nf90_close(fID),"close",filename2) -write(*,*) '"'//filename2//'" processed!' +write(*,*) '"> '//filename2//'" processed!' END SUBROUTINE read_data_PCM Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_mod.F90 (revision 3602) +++ trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_mod.F90 (revision 3603) @@ -40,5 +40,5 @@ real :: coef, avg -write(*,*) "Update of the CO2 tendency from the current pressure" +write(*,*) "> Updating the CO2 ice tendency for the new pressure" ! Evolution of the water ice for each physical point Index: trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90 (revision 3602) +++ trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90 (revision 3603) @@ -824,6 +824,6 @@ ierr=NF_INQ_DIMID(nid,"time",id(4)) ! Tell the world about it - write(*,*) "=====================" - write(*,*) "writediagsoilpem: creating variable "//trim(name) + write(*,*) "===========================" + write(*,*) "DIAGSOILPEM: creating variable "//trim(name) call def_var(nid,name,title,units,4,id,varid,ierr) endif ! of if (ierr.ne.NF_NOERR) @@ -908,6 +908,6 @@ ierr=NF_INQ_DIMID(nid,"time",id(3)) ! Tell the world about it - write(*,*) "=====================" - write(*,*) "writediagsoilpem: creating variable "//trim(name) + write(*,*) "===========================" + write(*,*) "DIAGSOILPEM: creating variable "//trim(name) call def_var(nid,name,title,units,3,id,varid,ierr) endif ! of if (ierr.ne.NF_NOERR) @@ -959,6 +959,6 @@ ierr=NF_INQ_DIMID(nid,"time",id(1)) ! Tell the world about it - write(*,*) "=====================" - write(*,*) "writediagsoilpem: creating variable "//trim(name) + write(*,*) "===========================" + write(*,*) "DIAGSOILPEM: creating variable "//trim(name) call def_var(nid,name,title,units,1,id,varid,ierr) endif ! of if (ierr.ne.NF_NOERR)