Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3142) +++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3143) @@ -145,2 +145,7 @@ == 22/11/2023 == JBC Update of files in the deftank: addition of variables in the xml definition files, inclusion of "callphys.def" into "run_PEM.def" and minor typo in "launch_pem.sh". + +== 29/11/2023 == JBC +'Watercap' has been removed from the water ice evolution since 'water_reservoir' does already the job + Some cleanings to simplify the code. + +/!\ Commit for the PEM management of h2o ice before a rework of ice management in the PEM! Index: /trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90 (revision 3142) +++ /trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90 (revision 3143) @@ -1,69 +1,76 @@ -module comsoil_h_PEM +MODULE comsoil_h_PEM implicit none - integer, parameter :: nsoilmx_PEM = 68 ! number of layers in the PEM - real,save,allocatable,dimension(:) :: layer_PEM ! soil layer depths [m] - real,save,allocatable,dimension(:) :: mlayer_PEM ! soil mid-layer depths [m] - real,save,allocatable,dimension(:,:,:) :: TI_PEM ! soil thermal inertia [SI] - real,save,allocatable,dimension(:,:) :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI] - ! variables (FC: built in firstcall in soil.F) - REAL,SAVE,ALLOCATABLE :: tsoil_PEM(:,:,:) ! sub-surface temperatures [K] - real,save,allocatable :: mthermdiff_PEM(:,:) ! (FC) mid-layer thermal diffusivity [SI] - real,save,allocatable :: thermdiff_PEM(:,:) ! (FC) inter-layer thermal diffusivity [SI] - real,save,allocatable :: coefq_PEM(:) ! (FC) q_{k+1/2} coefficients [SI] - real,save,allocatable :: coefd_PEM(:,:) ! (FC) d_k coefficients [SI] - real,save,allocatable :: alph_PEM(:,:) ! (FC) alpha_k coefficients [SI] - real,save,allocatable :: beta_PEM(:,:) ! beta_k coefficients [SI] - real,save :: mu_PEM ! mu coefficient [SI] - real,save :: fluxgeo ! Geothermal flux [W/m^2] - real,save :: depth_breccia ! Depth at which we have breccia [m] - real,save :: depth_bedrock ! Depth at which we have bedrock [m] - integer,save :: index_breccia ! last index of the depth grid before having breccia - integer,save :: index_bedrock ! last index of the depth grid before having bedrock - LOGICAL soil_pem ! True by default, to run with the subsurface physic. Read in pem.def - real,save,allocatable,dimension(:) :: water_reservoir ! Large reserve of water [kg/m^2] - real,save :: water_reservoir_nom ! Nominal value for the large reserve of water [kg/m^2] - logical, save :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure +integer, parameter :: nsoilmx_PEM = 68 ! number of layers in the PEM +real, allocatable, dimension(:), save :: layer_PEM ! soil layer depths [m] +real, allocatable, dimension(:), save :: mlayer_PEM ! soil mid-layer depths [m] +real, allocatable, dimension(:,:,:), save :: TI_PEM ! soil thermal inertia [SI] +real, allocatable, dimension(:,:), save :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI] +! Variables (FC: built in firstcall in soil.F) +real, allocatable, dimension(:,:,:), save :: tsoil_PEM ! sub-surface temperatures [K] +real, allocatable, dimension(:,:), save :: mthermdiff_PEM ! (FC) mid-layer thermal diffusivity [SI] +real, allocatable, dimension(:,:), save :: thermdiff_PEM ! (FC) inter-layer thermal diffusivity [SI] +real, allocatable, dimension(:), save :: coefq_PEM ! (FC) q_{k+1/2} coefficients [SI] +real, allocatable, dimension(:,:), save :: coefd_PEM ! (FC) d_k coefficients [SI] +real, allocatable, dimension(:,:), save :: alph_PEM ! (FC) alpha_k coefficients [SI] +real, allocatable, dimension(:,:), save :: beta_PEM ! beta_k coefficients [SI] +real, save :: mu_PEM ! mu coefficient [SI] +real, save :: fluxgeo ! Geothermal flux [W/m^2] +real, save :: depth_breccia ! Depth at which we have breccia [m] +real, save :: depth_bedrock ! Depth at which we have bedrock [m] +integer, save :: index_breccia ! last index of the depth grid before having breccia +integer, save :: index_bedrock ! last index of the depth grid before having bedrock +logical :: soil_pem ! True by default, to run with the subsurface physic. Read in pem.def +real, allocatable, dimension(:), save :: water_reservoir ! Large reserve of water [kg/m^2] +real, save :: water_reservoir_nom ! Nominal value for the large reserve of water [kg/m^2] +logical, save :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure + +!======================================================================= contains +!======================================================================= - subroutine ini_comsoil_h_PEM(ngrid,nslope) - - implicit none - integer,intent(in) :: ngrid ! number of atmospheric columns - integer,intent(in) :: nslope ! number of slope within a mesh +SUBROUTINE ini_comsoil_h_PEM(ngrid,nslope) - allocate(layer_PEM(nsoilmx_PEM)) - allocate(mlayer_PEM(0:nsoilmx_PEM-1)) - allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope)) - allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope)) - allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1)) - allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1)) - allocate(coefq_PEM(0:nsoilmx_PEM-1)) - allocate(coefd_PEM(ngrid,nsoilmx_PEM-1)) - allocate(alph_PEM(ngrid,nsoilmx_PEM-1)) - allocate(beta_PEM(ngrid,nsoilmx_PEM-1)) - allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) - allocate(water_reservoir(ngrid)) - end subroutine ini_comsoil_h_PEM +implicit none +integer, intent(in) :: ngrid ! number of atmospheric columns +integer, intent(in) :: nslope ! number of slope within a mesh - subroutine end_comsoil_h_PEM +allocate(layer_PEM(nsoilmx_PEM)) +allocate(mlayer_PEM(0:nsoilmx_PEM-1)) +allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope)) +allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope)) +allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1)) +allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1)) +allocate(coefq_PEM(0:nsoilmx_PEM-1)) +allocate(coefd_PEM(ngrid,nsoilmx_PEM-1)) +allocate(alph_PEM(ngrid,nsoilmx_PEM-1)) +allocate(beta_PEM(ngrid,nsoilmx_PEM-1)) +allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) +allocate(water_reservoir(ngrid)) - implicit none +END SUBROUTINE ini_comsoil_h_PEM - if (allocated(layer_PEM)) deallocate(layer_PEM) - if (allocated(mlayer_PEM)) deallocate(mlayer_PEM) - if (allocated(TI_PEM)) deallocate(TI_PEM) - if (allocated(tsoil_PEM)) deallocate(tsoil_PEM) - if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM) - if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM) - if (allocated(coefq_PEM)) deallocate(coefq_PEM) - if (allocated(coefd_PEM)) deallocate(coefd_PEM) - if (allocated(alph_PEM)) deallocate(alph_PEM) - if (allocated(beta_PEM)) deallocate(beta_PEM) - if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM) - if (allocated(water_reservoir)) deallocate(water_reservoir) - end subroutine end_comsoil_h_PEM +!======================================================================= -end module comsoil_h_PEM +SUBROUTINE end_comsoil_h_PEM + +implicit none + +if (allocated(layer_PEM)) deallocate(layer_PEM) +if (allocated(mlayer_PEM)) deallocate(mlayer_PEM) +if (allocated(TI_PEM)) deallocate(TI_PEM) +if (allocated(tsoil_PEM)) deallocate(tsoil_PEM) +if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM) +if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM) +if (allocated(coefq_PEM)) deallocate(coefq_PEM) +if (allocated(coefd_PEM)) deallocate(coefd_PEM) +if (allocated(alph_PEM)) deallocate(alph_PEM) +if (allocated(beta_PEM)) deallocate(beta_PEM) +if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM) +if (allocated(water_reservoir)) deallocate(water_reservoir) + +END SUBROUTINE end_comsoil_h_PEM + +END MODULE comsoil_h_PEM Index: /trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90 (revision 3142) +++ /trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90 (revision 3143) @@ -38,7 +38,6 @@ ! Conversion H2O/CO2 frost to perennial frost and vice versa -real, parameter :: threshold_water_frost2perennial = 1000. !~ 1 m -real, parameter :: threshold_co2_frost2perennial = 16000. !~ 10 m +real, parameter :: threshold_h2o_frost2perennial = 1000. !~ 1 m +real, parameter :: threshold_co2_frost2perennial = 16000. !~ 10 m END MODULE constants_marspem_mod - Index: /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90 (revision 3142) +++ /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90 (revision 3143) @@ -3,5 +3,7 @@ implicit none +! ---------------------------------------------------------------------- contains +! ---------------------------------------------------------------------- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -44,21 +46,28 @@ STOPPING = .false. -! Computation of the present surface of water ice sublimating +! Computation of the present surface of h2o ice present_surf = 0. do i = 1,ngrid do islope = 1,nslope - !if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope) - if (initial_h2o_ice(i,islope) > 0.5) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope) + if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope) enddo enddo ! Check of the criterion -if (present_surf < ini_surf*(1. - water_ice_criterion) .or. present_surf > ini_surf*(1. + water_ice_criterion)) then +if (present_surf < ini_surf*(1. - water_ice_criterion)) then STOPPING = .true. write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold" + write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", present_surf < ini_surf*(1. - water_ice_criterion) write(*,*) "Current surface of water ice sublimating =", present_surf write(*,*) "Initial surface of water ice sublimating =", ini_surf write(*,*) "Percentage of change accepted =", water_ice_criterion*100 - write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", (present_surf < ini_surf*(1. - water_ice_criterion)) +endif +if (present_surf > ini_surf*(1. + water_ice_criterion)) then + STOPPING = .true. + write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold" + write(*,*) "present_surf > ini_surf*(1. + water_ice_criterion)", present_surf > ini_surf*(1. + water_ice_criterion) + write(*,*) "Current surface of water ice sublimating =", present_surf + write(*,*) "Initial surface of water ice sublimating =", ini_surf + write(*,*) "Percentage of change accepted =", water_ice_criterion*100 endif @@ -105,5 +114,5 @@ STOPPING_ps = .false. -! Computation of the actual surface +! Computation of the present surface of co2 ice present_surf = 0. do i = 1,ngrid @@ -114,22 +123,38 @@ ! Check of the criterion -if (present_surf < ini_surf*(1. - co2_ice_criterion) .or. present_surf > ini_surf*(1. + co2_ice_criterion)) then +if (present_surf < ini_surf*(1. - co2_ice_criterion)) then STOPPING_ice = .true. - write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold" + write(*,*) "Reason of stopping: the surface of co2 ice sublimating reached the threshold" + write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", present_surf < ini_surf*(1. - co2_ice_criterion) write(*,*) "Current surface of co2 ice sublimating =", present_surf write(*,*) "Initial surface of co2 ice sublimating =", ini_surf write(*,*) "Percentage of change accepted =", co2_ice_criterion*100. - write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", (present_surf < ini_surf*(1. - co2_ice_criterion)) +endif +if (present_surf > ini_surf*(1. + co2_ice_criterion)) then + STOPPING_ice = .true. + write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold" + write(*,*) "present_surf > ini_surf*(1. + co2_ice_criterion)", present_surf > ini_surf*(1. + co2_ice_criterion) + write(*,*) "Current surface of co2 ice sublimating =", present_surf + write(*,*) "Initial surface of co2 ice sublimating =", ini_surf + write(*,*) "Percentage of change accepted =", co2_ice_criterion*100. endif if (abs(ini_surf) < 1.e-5) STOPPING_ice = .false. -if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion) .or. global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then +if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)) then STOPPING_ps = .true. write(*,*) "Reason of stopping: the global pressure reach the threshold" + write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion) write(*,*) "Current global pressure =", global_ave_press_new write(*,*) "PCM global pressure =", global_ave_press_GCM write(*,*) "Percentage of change accepted =", ps_criterion*100. - write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)) +endif +if (global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then + STOPPING_ps = .true. + write(*,*) "Reason of stopping: the global pressure reach the threshold" + write(*,*) "global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)", global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion) + write(*,*) "Current global pressure =", global_ave_press_new + write(*,*) "PCM global pressure =", global_ave_press_GCM + write(*,*) "Percentage of change accepted =", ps_criterion*100. endif Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3142) +++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3143) @@ -1,10 +1,10 @@ !------------------------ ! I Initialization -! I_a READ run.def -! I_b READ of start_evol.nc and starfi_evol.nc +! I_a Read run.def +! I_b Read of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation -! I_d READ PCM data and convert to the physical grid +! I_d Read PCM data and convert to the physical grid ! I_e Initialization of the PEM variable and soil -! I_f Compute tendencies & Save initial situation +! I_f Compute tendencies ! I_g Save initial PCM situation ! I_h Read the startpem.nc @@ -42,5 +42,5 @@ use criterion_pem_stop_mod, only: criterion_waterice_stop, criterion_co2_stop use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, & - m_noco2, threshold_water_frost2perennial, threshold_co2_frost2perennial + m_noco2, threshold_h2o_frost2perennial, threshold_co2_frost2perennial use evol_co2_ice_s_mod, only: evol_co2_ice_s use evol_h2o_ice_s_mod, only: evol_h2o_ice_s @@ -86,4 +86,5 @@ use paleoclimate_mod, only: albedo_perennialco2 use comcstfi_h, only: pi, rad, g, cpp, mugaz, r + use surfini_mod, only: surfini #else use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names @@ -131,5 +132,5 @@ ! Variables to read start.nc -character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM +character(*), parameter :: start_name = "start_evol.nc" ! Name of the file used to initialize the PEM ! Dynamic variables @@ -138,7 +139,7 @@ real, dimension(ip1jmp1,llm) :: teta ! temperature potentielle real, dimension(:,:,:), allocatable :: q ! champs advectes -real, dimension(ip1jmp1) :: ps ! pression au sol +real, dimension(ip1jmp1) :: ps ! pression au sol real, dimension(:), allocatable :: ps_start_GCM ! (ngrid) pression au sol -real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées +real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées real, dimension(ip1jmp1,llm) :: masse ! masse d'air real, dimension(ip1jmp1) :: phis ! geopotentiel au sol @@ -146,15 +147,15 @@ ! Variables to read starfi.nc -character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" ! Name of the file used for initialsing the PEM -character(2) :: str2 -integer :: ncid, varid, status ! Variable for handling opening of files -integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files -integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files -integer :: apvarid, bpvarid ! Variable ID for Netcdf files +character(*), parameter :: startfi_name = "startfi_evol.nc" ! Name of the file used to initialize the PEM +character(2) :: str2 +integer :: ncid, varid, status ! Variable for handling opening of files +integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files +integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files +integer :: apvarid, bpvarid ! Variable ID for Netcdf files ! Variables to read starfi.nc and write restartfi.nc -real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAME and written in restartfi -real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi -real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi +real, dimension(:), allocatable :: longitude ! Longitude read in startfi_name and written in restartfi +real, dimension(:), allocatable :: latitude ! Latitude read in startfi_name and written in restartfi +real, dimension(:), allocatable :: cell_area ! Cell_area read in startfi_name and written in restartfi real :: Total_surface ! Total surface of the planet @@ -187,6 +188,4 @@ real, dimension(:,:), allocatable :: co2frost_ini ! Initial amount of co2 frost (at the start of the PEM run) real, dimension(:,:), allocatable :: perennial_co2ice_ini ! Initial amoun of perennial co2 ice (at the start of the PEM run) - - ! Variables for slopes @@ -243,9 +242,9 @@ real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice real :: albedo_h2o_frost ! albedo of h2o frost - real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic - real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic - real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic - real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic - real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic + real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic + real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic + real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic + real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic + real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case @@ -259,7 +258,7 @@ #ifdef CPP_1D - integer :: nsplit_phys - integer, parameter :: jjm_value = jjm - 1 - integer :: ierr + integer :: nsplit_phys + integer, parameter :: jjm_value = jjm - 1 + integer :: ierr ! Dummy variables to use the subroutine 'init_testphys1d' @@ -272,7 +271,7 @@ real, dimension(nlayer + 1) :: plev #else - integer, parameter :: jjm_value = jjm - real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart - real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart + integer, parameter :: jjm_value = jjm + real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart + real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart #endif @@ -302,5 +301,5 @@ !------------------------ ! I Initialization -! I_a READ run.def +! I_a Read run.def !------------------------ #ifndef CPP_1D @@ -315,5 +314,5 @@ allocate(longitude(1),latitude(1),cell_area(1)) - therestart1D = .false. + therestart1D = .false. ! Default value inquire(file = 'start1D_evol.txt',exist = therestart1D) if (.not. therestart1D) then @@ -321,5 +320,5 @@ error stop 'Initialization cannot be done for the 1D PEM.' endif - therestartfi = .false. + therestartfi = .false. ! Default value inquire(file = 'startfi_evol.nc',exist = therestartfi) if (.not. therestartfi) then @@ -339,10 +338,10 @@ !------------------------ ! I Initialization -! I_b READ of start_evol.nc and starfi_evol.nc -!------------------------ -! I_b.1 READ start_evol.nc +! I_b Read of start_evol.nc and starfi_evol.nc +!------------------------ +! I_b.1 Read start_evol.nc allocate(ps_start_GCM(ngrid)) #ifndef CPP_1D - call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0) + call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0) call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM) @@ -353,5 +352,5 @@ allocate(ap(nlayer + 1)) allocate(bp(nlayer + 1)) - status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid) + status = nf90_open(start_name,NF90_NOWRITE,ncid) status = nf90_inq_varid(ncid,"ap",apvarid) status = nf90_get_var(ncid,apvarid,ap) @@ -368,5 +367,5 @@ ! In the PCM, these values are given to the physic by the dynamic. ! Here we simply read them in the startfi_evol.nc file -status = nf90_open(FILE_NAME, NF90_NOWRITE, ncid) +status = nf90_open(startfi_name, NF90_NOWRITE, ncid) status = nf90_inq_varid(ncid,"longitude",lonvarid) @@ -384,8 +383,8 @@ status = nf90_close(ncid) -! I_b.2 READ startfi_evol.nc +! I_b.2 Read startfi_evol.nc ! First we read the initial state (starfi.nc) #ifndef CPP_STD - call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, & + 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, & watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist) @@ -394,5 +393,5 @@ where (qsurf < 0.) qsurf = 0. - call surfini(ngrid,qsurf) + call surfini(ngrid,nslope,qsurf) #else call phys_state_var_init(nq) @@ -414,5 +413,5 @@ allocate(albedo(ngrid,2,1)) allocate(inertiesoil(ngrid,nsoilmx,1)) - call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, & + call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, & tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, & qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, & @@ -449,5 +448,5 @@ if (noms(nnq) == "h2o_vap") then igcm_h2o_vap = nnq - mmol(igcm_h2o_vap)=18. + mmol(igcm_h2o_vap) = 18. endif if (noms(nnq) == "co2") igcm_co2 = nnq @@ -480,5 +479,5 @@ !------------------------ ! I Initialization -! I_d READ PCM data and convert to the physical grid +! I_d Read PCM data and convert to the physical grid !------------------------ ! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value @@ -552,5 +551,5 @@ !------------------------ ! I Initialization -! I_f Compute tendencies & Save initial situation +! I_f Compute tendencies !------------------------ allocate(tendencies_co2_ice(ngrid,nslope)) @@ -603,5 +602,5 @@ write(*,*) "Total initial surface of co2 ice sublimating (slope) =", ini_surf_co2 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o -write(*,*) "Total surface of the planet=", Total_surface +write(*,*) "Total surface of the planet =", Total_surface allocate(zplev_gcm(ngrid,nlayer + 1)) @@ -617,5 +616,5 @@ global_ave_press_GCM = global_ave_press_old global_ave_press_new = global_ave_press_old -write(*,*) "Initial global average pressure=", global_ave_press_GCM +write(*,*) "Initial global average pressure =", global_ave_press_GCM !------------------------ @@ -638,5 +637,5 @@ do ig = 1,ngrid do islope = 1,nslope - qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) + qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) + perennial_co2ice(ig,islope) enddo @@ -657,6 +656,6 @@ enddo - write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded - write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded + write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded + write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded endif ! adsorption deallocate(tsurf_ave_yr1) @@ -735,6 +734,9 @@ q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ & (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) - if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30 - if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1. + if (q_h2o_PEM_phys(ig,t) < 0) then + q_h2o_PEM_phys(ig,t) = 1.e-30 + else if (q_h2o_PEM_phys(ig,t) > 1) then + q_h2o_PEM_phys(ig,t) = 1. + endif enddo enddo @@ -749,5 +751,5 @@ if (q_co2_PEM_phys(ig,t) < 0) then q_co2_PEM_phys(ig,t) = 1.e-30 - elseif (q_co2_PEM_phys(ig,t) > 1) then + else if (q_co2_PEM_phys(ig,t) > 1) then q_co2_PEM_phys(ig,t) = 1. endif @@ -898,5 +900,5 @@ call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM) -! II_d.5 Update the mass of the regolith adsorbded +! II_d.5 Update the mass of the regolith adsorbed if (adsorption_pem) then call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, & @@ -985,5 +987,5 @@ exit else - write(*,*) "We continue!" + write(*,*) "The PEM can continue!" write(*,*) "Number of iterations of the PEM: year_iter =", year_iter write(*,*) "Number of simulated Martian years: i_myear =", i_myear @@ -1007,9 +1009,11 @@ watercap_sum = 0. do islope = 1,nslope - if (qsurf(ig,igcm_h2o_ice,islope) > (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy - qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) ! put back ancien frost + ! water_reservoir and watercap have not changed since PCM call: here we check if we have accumulate frost or not. + if (qsurf(ig,igcm_h2o_ice,islope) > water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) then + ! 1st case: more ice than initialy. Put back ancien frost + qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) else -! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap - watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) + ! 2nd case: ice sublimated. Then let's put qsurf = 0 and add the difference in watercap + watercap(ig,islope) = qsurf(ig,igcm_h2o_ice,islope) qsurf(ig,igcm_h2o_ice,islope) = 0. endif @@ -1027,13 +1031,13 @@ enddo if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. - if (water_sum > threshold_water_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1 + if (water_sum > threshold_h2o_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1 watercaptag(ig) = .true. - water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost + water_reservoir(ig) = water_reservoir(ig) + threshold_h2o_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost do islope = 1,nslope - qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.) + qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_h2o_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.) enddo endif else ! let's check that the infinite source of water disapear - if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perennial) then + if ((water_sum + water_reservoir(ig)) < threshold_h2o_frost2perennial) then watercaptag(ig) = .false. do islope = 1,nslope Index: /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3142) +++ /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3143) @@ -104,5 +104,5 @@ !0. Check if the start_PEM exist. -inquire(file=filename,exist = startpem_file) +inquire(file = filename,exist = startpem_file) write(*,*)'Is start PEM?',startpem_file @@ -312,15 +312,10 @@ !. 5 water reservoir #ifndef CPP_STD + water_reservoir = 0. call get_field("water_reservoir",water_reservoir,found) - if(.not.found) then - write(*,*)'Pemetat0: failed loading ' - write(*,*)'will reconstruct the values from watercaptag' - do ig=1,ngrid - if(watercaptag(ig)) then - water_reservoir(ig)=water_reservoir_nom - else - water_reservoir(ig)=0. - endif - enddo + if (.not. found) then + write(*,*)'Pemetat0: failed loading ' + write(*,*)'will reconstruct the values from watercaptag' + where (watercaptag) water_reservoir = water_reservoir_nom endif #endif @@ -481,11 +476,6 @@ !. e) water reservoir #ifndef CPP_STD - do ig = 1,ngrid - if (watercaptag(ig)) then - water_reservoir(ig) = water_reservoir_nom - else - water_reservoir(ig) = 0. - endif - enddo + water_reservoir = 0. + where (watercaptag) water_reservoir = water_reservoir_nom #endif Index: /trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 (revision 3142) +++ /trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 (revision 3143) @@ -104,4 +104,7 @@ write(*,*) "Data for h2o_ice_s downloaded!" + call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:)) + write(*,*) "Data for tsurf downloaded!" + #ifndef CPP_STD call get_var3("watercap",watercap(:,:,1,:)) @@ -110,10 +113,5 @@ call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:)) write(*,*) "Data for perennial_co2ice downloaded!" -#endif - - call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:)) - write(*,*) "Data for tsurf downloaded!" - -#ifndef CPP_STD + if (soil_pem) then call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:)) @@ -140,4 +138,10 @@ write(*,*) "Data for h2o_ice_s downloaded!" + do islope = 1,nslope + write(num,'(i2.2)') islope + call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:)) + enddo + write(*,*) "Data for tsurf downloaded!" + #ifndef CPP_STD do islope = 1,nslope @@ -152,13 +156,5 @@ enddo write(*,*) "Data for perennial_co2ice downloaded!" -#endif - - do islope = 1,nslope - write(num,'(i2.2)') islope - call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:)) - enddo - write(*,*) "Data for tsurf downloaded!" - -#ifndef CPP_STD + if (soil_pem) then do islope = 1,nslope Index: /trunk/LMDZ.MARS/deftank/pem/concat_diagpem.sh =================================================================== --- /trunk/LMDZ.MARS/deftank/pem/concat_diagpem.sh (revision 3142) +++ /trunk/LMDZ.MARS/deftank/pem/concat_diagpem.sh (revision 3143) @@ -58,3 +58,3 @@ done -echo "Concatenation of \"$directory/diagpem*.nc\" files into \"$directory/$output_file\" is complete!" +echo "Concatenation of \"$directory/diagpem*.nc\" files into \"$output_file\" is complete!"