Changeset 3603

Timestamp:

Jan 28, 2025, 5:16:15 PM (37 hours ago)

Author:

jbclement

Message:

PEM:

• 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.

JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• evol_ice_mod.F90 (modified) (2 diffs)
• glaciers_mod.F90 (modified) (2 diffs)
• pem.F90 (modified) (35 diffs)
• read_data_PCM_mod.F90 (modified) (2 diffs)
• recomp_tend_co2_mod.F90 (modified) (1 diff)
• writediagpem.F90 (modified) (3 diffs)

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

@@ -562 +562 @@
 - 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.

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

@@ -34 +34 @@
 !=======================================================================
 ! Evolution of CO2 ice for each physical point
-write(*,*) 'Evolution of co2 ice'
+write(*,*) '> Evolution of CO2 ice'
 
 co2_ice_old = co2_ice
@@ -88 +88 @@
 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

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

@@ -53 +53 @@
 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 +93 @@
 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)

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

@@ -139 +139 @@
 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 +286 @@
 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 +380 @@
 !    I_a Read the "run.def"
 !------------------------
+write(*,*) '> Reading "run.def" (PEM)'
 #ifndef CPP_1D
     dtphys = daysec/48. ! Dummy value (overwritten in phyetat0)
@@ -379 +389 @@
     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 +405 @@
     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 +420 @@
 !------------------------
 ! 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 +443 @@
     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 +463 @@
 ! 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 +524 @@
 #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 +543 @@
 enddo
 write(*,*) 'Flat slope for islope = ',iflat
-write(*,*) 'corresponding criterium = ',def_slope_mean(iflat)
+write(*,*) 'Corresponding criterium = ',def_slope_mean(iflat)
 
 !------------------------
@@ -597 +611 @@
 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 +618 @@
 !    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 +657 @@
     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 +721 @@
 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 +734 @@
         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 +745 @@
         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 +757 @@
 
 ! 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 +818 @@
 !------------------------
 ! II_d.1 Update Tsurf
-    write(*,*) "Updating the new Tsurf"
+    write(*,*) "> Updating surface temperature"
     do ig = 1,ngrid
         do islope = 1,nslope
@@ -811 +824 @@
             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 +851 @@
     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 +877 @@
         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 +929 @@
                 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 +978 @@
 !    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 +1011 @@
     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 +1037 @@
 !------------------------
 ! III_a.1 Ice update for start file
+write(*,*) '> Reconstructing perennial ice and frost for the PCM'
 watercap = 0.
 perennial_co2ice = co2_ice
@@ -1052 +1066 @@
 
 ! 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 +1072 @@
 ! 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 +1086 @@
 
 ! 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 +1092 @@
 
 ! 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 +1136 @@
 
 ! 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 +1143 @@
         icap = 1 ! Northern hemisphere
     endif
-    do islope = 1,ngrid
+    do islope = 1,nslope
         ! Bare ground
         albedo(ig,:,islope) = albedodat(ig)
@@ -1154 +1173 @@
 
 ! 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 +1184 @@
 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 +1228 @@
                   tsea_ice,sea_ice)
 #endif
-write(*,*) "restartfi.nc has been written."
 
 !------------------------
@@ -1213 +1233 @@
 !    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

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

@@ -181 +181 @@
 ! Close the NetCDF file
 call error_msg(nf90_close(fID),"close",filename1)
-write(*,*) '"'//filename1//'" processed!'
+write(*,*) '> "'//filename1//'" processed!'
 
 !------------------------------- Year 2 --------------------------------
@@ -445 +445 @@
 ! Close the NetCDF file
 call error_msg(nf90_close(fID),"close",filename2)
-write(*,*) '"'//filename2//'" processed!'
+write(*,*) '"> '//filename2//'" processed!'
 
 END SUBROUTINE read_data_PCM

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

@@ -40 +40 @@
 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

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

@@ -824 +824 @@
     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 +908 @@
     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 +959 @@
     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)