Changeset 3130 for trunk/LMDZ.COMMON/libf
- Timestamp:
- Nov 20, 2023, 1:25:31 PM (13 months ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90
r3082 r3130 156 156 #ifndef CPP_STD 157 157 158 !0.1 Look at peren ial ice158 !0.1 Look at perennial ice 159 159 do ig = 1,ngrid 160 160 do islope = 1,nslope … … 317 317 #ifndef CPP_STD 318 318 319 !0.1 Look at peren ial ice319 !0.1 Look at perennial ice 320 320 do ig = 1,ngrid 321 321 do islope = 1,nslope -
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3123 r3130 137 137 Adapting the PEM soil grid to the one of the PCM 138 138 Minor corrections when reading/initializing soil temperature, subsurface water ice 139 140 == 20/11/2023 == JBC 141 The perennial co2 ice is now taken into account with co2 frost (qsurf) to compute the tendency and to make the update + Rework of how co2 frost is converted to perennial co2 ice at the end of the PEM run + Correction of the value of 'threshold_co2_frost2perennial' to correspond to 10 m + Perennial co2 ice is now handled outside 'paleoclimate' in "phyetat0_mod.F90" of the Mars PCM + Some cleanings. 142 143 /!\ Commit for the PEM management of co2 ice before a rework of ice management in the PEM! -
trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope_mod.F90
r3076 r3130 25 25 26 26 ! OUTPUT 27 real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of peren ial ice27 real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of perennial ice 28 28 !======================================================================= 29 29 -
trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90
r3076 r3130 37 37 real, parameter :: Lco2 = 5.71e5 ! Pilorget and Forget 2016 38 38 39 ! Conversion H2O/CO2 frost to peren ial frost and vice versa40 real, parameter :: threshold_water_frost2peren ial = 1000. !~ 1m41 real, parameter :: threshold_co2_frost2peren ial = 3200. !~ 2m39 ! Conversion H2O/CO2 frost to perennial frost and vice versa 40 real, parameter :: threshold_water_frost2perennial = 1000. !~ 1 m 41 real, parameter :: threshold_co2_frost2perennial = 16000. !~ 10 m 42 42 43 43 END MODULE constants_marspem_mod -
trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90
r3050 r3130 1 module criterion_pem_stop_mod 2 implicit none 1 MODULE criterion_pem_stop_mod 3 2 4 contains 3 implicit none 5 4 5 contains 6 6 7 7 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 8 8 !!! 9 !!! Purpose: Criterions to check if the PEM needs to call the GCM !!!9 !!! Purpose: Criterions to check if the PEM needs to call the PCM 10 10 !!! Author: RV & LL, 02/2023 11 11 !!! 12 12 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 13 13 14 15 14 SUBROUTINE criterion_waterice_stop(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice) 16 15 17 18 use comslope_mod, only: subslope_dist,nslope16 use time_evol_mod, only: water_ice_criterion 17 use comslope_mod, only: subslope_dist, nslope 19 18 20 IMPLICIT NONE 19 implicit none 21 20 22 21 !======================================================================= … … 28 27 ! arguments: 29 28 ! ---------- 30 31 29 ! INPUT 32 INTEGER, intent(in) :: ngrid ! # of grid physical grid points 33 REAL, intent(in) :: cell_area(ngrid) ! physical point field : Area of the cells 34 REAL, intent(in) :: qsurf(ngrid,nslope) ! physical point field : Actual density of water ice 35 REAL, intent(in) :: ini_surf ! Initial surface of h2o ice that was sublimating 36 REAL, intent(in) :: initial_h2o_ice(ngrid,nslope) ! Grid point that initialy were covered by h2o_ice 37 30 integer, intent(in) :: ngrid ! # of physical grid points 31 real, dimension(ngrid), intent(in) :: cell_area ! Area of the cells 32 real, dimension(ngrid,nslope), intent(in) :: qsurf ! Actual density of water ice 33 real, intent(in) :: ini_surf ! Initial surface of h2o ice that was sublimating 34 real, dimension(ngrid,nslope), intent(in) :: initial_h2o_ice ! Grid point that initialy were covered by h2o_ice 38 35 ! OUTPUT 39 LOGICAL, intent(out) :: STOPPING ! Logical : is the criterion reached? 40 36 logical, intent(out) :: STOPPING ! Is the criterion reached? 41 37 ! local: 42 ! ----- 43 INTEGER :: i,islope! Loop44 REAL :: present_surf! Initial/Actual surface of water ice38 ! ------ 39 integer :: i, islope ! Loop 40 real :: present_surf ! Initial/Actual surface of water ice 45 41 46 42 !======================================================================= 43 ! Initialisation to false 44 STOPPING = .false. 47 45 48 ! initialisation to false 49 STOPPING=.FALSE. 46 ! Computation of the present surface of water ice sublimating 47 present_surf = 0. 48 do i = 1,ngrid 49 do islope = 1,nslope 50 !if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope) 51 if (initial_h2o_ice(i,islope) > 0.5) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope) 52 enddo 53 enddo 50 54 51 ! computation of the present surface of water ice sublimating 52 present_surf=0. 53 do i=1,ngrid 54 do islope=1, nslope 55 if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then 56 present_surf=present_surf+cell_area(i)*subslope_dist(i,islope) 57 endif 58 enddo 59 enddo 60 61 ! check of the criterion 62 if(present_surf.LT.ini_surf*(1-water_ice_criterion) .OR. & 63 present_surf.GT.ini_surf*(1+water_ice_criterion)) then 64 STOPPING=.TRUE. 65 write(*,*) "Reason of stopping : The surface of water ice sublimating reach the threshold:" 66 write(*,*) "Current surface of water ice sublimating=", present_surf 67 write(*,*) "Initial surface of water ice sublimating=", ini_surf 68 write(*,*) "Percentage of change accepted=", water_ice_criterion*100 69 write(*,*) "present_surf<ini_surf*(1-water_ice_criterion)", (present_surf.LT.ini_surf*(1-water_ice_criterion)) 70 endif 55 ! Check of the criterion 56 if (present_surf < ini_surf*(1. - water_ice_criterion) .or. present_surf > ini_surf*(1. + water_ice_criterion)) then 57 STOPPING = .true. 58 write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold" 59 write(*,*) "Current surface of water ice sublimating =", present_surf 60 write(*,*) "Initial surface of water ice sublimating =", ini_surf 61 write(*,*) "Percentage of change accepted =", water_ice_criterion*100 62 write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", (present_surf < ini_surf*(1. - water_ice_criterion)) 63 endif 71 64 72 if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then 73 STOPPING=.FALSE. 74 endif 65 if (ini_surf < 1.e-5 .and. ini_surf > -1.e-5) STOPPING = .false. 66 75 67 END SUBROUTINE criterion_waterice_stop 76 68 77 ! ---------------------------------------------------------------------- --------------------------69 ! ---------------------------------------------------------------------- 78 70 79 71 SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_co2_ice,global_ave_press_GCM,global_ave_press_new,nslope) 80 72 81 use time_evol_mod, only: co2_ice_criterion,ps_criterion82 73 use time_evol_mod, only: co2_ice_criterion, ps_criterion 74 use comslope_mod, only: subslope_dist 83 75 84 IMPLICIT NONE 76 implicit none 85 77 86 78 !======================================================================= 87 79 ! 88 ! Routine that checks if the c riterion to stop the PEM isreached80 ! Routine that checks if the co2 and pressure criteria to stop the PEM are reached 89 81 ! 90 82 !======================================================================= … … 92 84 ! arguments: 93 85 ! ---------- 94 95 86 ! INPUT 96 INTEGER, intent(in) :: ngrid,nslope ! # of grid physical grid points 97 REAL, intent(in) :: cell_area(ngrid) ! physical point field : Area of the cells 98 REAL, intent(in) :: qsurf(ngrid,nslope) ! physical point field : Actual density of water ice 99 REAL, intent(in) :: ini_surf ! Initial surface of co2 ice that was sublimating 100 REAL, intent(in) :: initial_co2_ice(ngrid,nslope) ! Grid point that initialy were covered by co2_ice 101 REAL, intent(in) :: global_ave_press_GCM ! Planet average pressure from the GCM start files 102 REAL, intent(in) :: global_ave_press_new ! Planet average pressure from the PEM computations 103 87 integer, intent(in) :: ngrid, nslope ! # of grid physical grid points 88 real, dimension(ngrid), intent(in) :: cell_area ! Area of the cells 89 real, dimension(ngrid,nslope), intent(in) :: qsurf ! Actual density of water ice 90 real, intent(in) :: ini_surf ! Initial surface of co2 ice that was sublimating 91 real, dimension(ngrid,nslope), intent(in) :: initial_co2_ice ! Grid point that initialy were covered by co2_ice 92 real, intent(in) :: global_ave_press_GCM ! Planet average pressure from the PCM start files 93 real, intent(in) :: global_ave_press_new ! Planet average pressure from the PEM computations 104 94 ! OUTPUT 105 LOGICAL, intent(out) :: STOPPING_ice ! Logical : is the criterion forice reached?106 LOGICAL, intent(out) :: STOPPING_ps ! Logical : is the criterion for pressure reached?95 logical, intent(out) :: STOPPING_ice ! Is the criterion for co2 ice reached? 96 logical, intent(out) :: STOPPING_ps ! Is the criterion for pressure reached? 107 97 ! local: 108 ! ----- 109 INTEGER :: i,islope! Loop110 REAL :: present_surf! Initial/Actual surface of water ice98 ! ------ 99 integer :: i, islope ! Loop 100 real :: present_surf ! Initial/Actual surface of water ice 111 101 112 102 !======================================================================= 103 ! Initialisation to false 104 STOPPING_ice = .false. 105 STOPPING_ps = .false. 113 106 114 ! initialisation to false 115 STOPPING_ice=.FALSE. 116 STOPPING_ps =.FALSE. 117 ! computation of the actual surface 118 present_surf=0. 119 do i=1,ngrid 120 do islope=1,nslope 121 if (initial_co2_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then 122 present_surf=present_surf+cell_area(i)*subslope_dist(i,islope) 123 endif 124 enddo 125 enddo 126 127 ! check of the criterion 128 if(present_surf.LT.ini_surf*(1-co2_ice_criterion) .OR. & 129 present_surf.GT.ini_surf*(1+co2_ice_criterion)) then 130 STOPPING_ice=.TRUE. 131 write(*,*) "Reason of stopping : The surface of co2 ice sublimating reach the threshold:" 132 write(*,*) "Current surface of co2 ice sublimating=", present_surf 133 write(*,*) "Initial surface of co2 ice sublimating=", ini_surf 134 write(*,*) "Percentage of change accepted=", co2_ice_criterion*100 135 write(*,*) "present_surf<ini_surf*(1-co2_ice_criterion)", (present_surf.LT.ini_surf*(1-co2_ice_criterion)) 136 endif 107 ! Computation of the actual surface 108 present_surf = 0. 109 do i = 1,ngrid 110 do islope = 1,nslope 111 if (initial_co2_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope) 112 enddo 113 enddo 137 114 138 if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then 139 STOPPING_ice=.FALSE. 140 endif 115 ! Check of the criterion 116 if (present_surf < ini_surf*(1. - co2_ice_criterion) .or. present_surf > ini_surf*(1. + co2_ice_criterion)) then 117 STOPPING_ice = .true. 118 write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold" 119 write(*,*) "Current surface of co2 ice sublimating =", present_surf 120 write(*,*) "Initial surface of co2 ice sublimating =", ini_surf 121 write(*,*) "Percentage of change accepted =", co2_ice_criterion*100. 122 write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", (present_surf < ini_surf*(1. - co2_ice_criterion)) 123 endif 141 124 142 if(global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion) .OR. & 143 global_ave_press_new.GT.global_ave_press_GCM*(1+ps_criterion)) then 144 STOPPING_ps=.TRUE. 145 write(*,*) "Reason of stopping : The global pressure reach the threshold:" 146 write(*,*) "Current global pressure=", global_ave_press_new 147 write(*,*) "GCM global pressure=", global_ave_press_GCM 148 write(*,*) "Percentage of change accepted=", ps_criterion*100 149 write(*,*) "global_ave_press_new<global_ave_press_GCM*(ps_criterion)", (global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion)) 150 endif 125 if (abs(ini_surf) < 1.e-5) STOPPING_ice = .false. 126 127 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 128 STOPPING_ps = .true. 129 write(*,*) "Reason of stopping: the global pressure reach the threshold" 130 write(*,*) "Current global pressure =", global_ave_press_new 131 write(*,*) "PCM global pressure =", global_ave_press_GCM 132 write(*,*) "Percentage of change accepted =", ps_criterion*100. 133 write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)) 134 endif 151 135 152 136 END SUBROUTINE criterion_co2_stop 153 137 154 155 138 END MODULE -
trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod.F90
r3039 r3130 1 1 MODULE evol_co2_ice_s_mod 2 2 3 IMPLICIT NONE 3 implicit none 4 4 5 CONTAINS 5 !======================================================================= 6 contains 7 !======================================================================= 6 8 7 SUBROUTINE evol_co2_ice_s(qsurf,tendencies_co2_ice_phys,& 8 iim_input,jjm_input,ngrid,cell_area,nslope) 9 SUBROUTINE evol_co2_ice_s(qsurf,tendencies_co2_ice_phys,iim_input,jjm_input,ngrid,cell_area,nslope) 9 10 10 11 use time_evol_mod, only: dt_pem 11 12 12 IMPLICIT NONE 13 implicit none 13 14 14 15 !======================================================================= … … 17 18 ! 18 19 !======================================================================= 19 20 20 ! arguments: 21 21 ! ---------- 22 23 22 ! INPUT 24 25 INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope ! # of grid points along longitude/latitude/ total 26 REAL, intent(in) :: cell_area(ngrid) 23 integer, intent(in) :: iim_input, jjm_input, ngrid, nslope ! # of grid points along longitude/latitude/ total 24 REAL, dimension(ngrid), intent(in) :: cell_area 27 25 28 26 ! OUTPUT 29 REAL, INTENT(INOUT) :: qsurf(ngrid,nslope) ! physical point field: Previous and actual density of water ice30 REAL, intent(inout) :: tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year27 real, dimension(ngrid,nslope), intent(inout) :: qsurf ! physical point field: Previous and actual density of water ice 28 real, dimension(ngrid,nslope), intent(inout) :: tendencies_co2_ice_phys ! physical point field: Evolution of perennial ice over one year 31 29 32 30 ! local: 33 31 ! ---- 34 35 INTEGER :: i,j,ig0,islope ! loop variable 36 32 integer :: i, j, ig0, islope ! loop variable 33 !======================================================================= 37 34 ! Evolution of the CO2 ice for each physical point 38 35 do i=1,ngrid 39 36 do islope=1,nslope 40 qsurf(i,islope)=qsurf(i,islope)+tendencies_co2_ice_phys(i,islope)*dt_pem41 if (qsurf(i,islope).lt.0) then42 qsurf(i,islope)=0.43 tendencies_co2_ice_phys(i,islope)=0.44 endif37 qsurf(i,islope) = qsurf(i,islope) + tendencies_co2_ice_phys(i,islope)*dt_pem 38 if (qsurf(i,islope) < 0) then 39 qsurf(i,islope) = 0. 40 tendencies_co2_ice_phys(i,islope) = 0. 41 endif 45 42 enddo 46 43 enddo 47 44 48 45 END SUBROUTINE evol_co2_ice_s -
trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod.F90
r3082 r3130 1 1 MODULE evol_h2o_ice_s_mod 2 2 3 IMPLICIT NONE 3 implicit none 4 4 5 CONTAINS 5 !======================================================================= 6 contains 7 !======================================================================= 6 8 7 9 SUBROUTINE evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf,tendencies_h2o_ice_phys,STOPPING) 8 10 9 use time_evol_mod, only: dt_pem 10 use comslope_mod, only: subslope_dist, def_slope_mean 11 use criterion_pem_stop_mod, only: criterion_waterice_stop 11 use time_evol_mod, only: dt_pem 12 use comslope_mod, only: subslope_dist, def_slope_mean 13 use criterion_pem_stop_mod, only: criterion_waterice_stop 14 12 15 #ifndef CPP_STD 13 16 use comcstfi_h, only: pi … … 16 19 #endif 17 20 18 IMPLICIT NONE 21 implicit none 19 22 20 23 !======================================================================= … … 23 26 ! 24 27 !======================================================================= 25 26 28 ! arguments: 27 29 ! ---------- 28 29 30 ! INPUT 30 31 INTEGER, intent(in) :: ngrid ! # of grid points along longitude/latitude grid; 32 INTEGER, intent(in) :: nslope ! # of subslope 33 REAL, intent(in) :: cell_area(ngrid) ! Area of each mesh grid (m^2) 34 REAL, intent(in) :: delta_h2o_adsorbded(ngrid) ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2) 35 REAL, intent(in) :: delta_h2o_icetablesublim(ngrid) ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2) 31 integer, intent(in) :: ngrid ! # of grid points along longitude/latitude grid; 32 integer, intent(in) :: nslope ! # of subslope 33 real, dimension(ngrid), intent(in) :: cell_area ! Area of each mesh grid (m^2) 34 real, dimension(ngrid), intent(in) :: delta_h2o_adsorbded ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2) 35 real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2) 36 36 37 37 ! OUTPUT 38 REAL, INTENT(INOUT) :: qsurf(ngrid,nslope) ! physical point field: Previous and actual density of water ice (kg/m^2)39 REAL, intent(inout) :: tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year (kg/m^2/year)40 LOGICAL, INTENT(INOUT) :: STOPPING! Stopping criterion38 real, dimension(ngrid,nslope), intent(inout) :: qsurf ! physical point field: Previous and actual density of water ice (kg/m^2) 39 real, dimension(ngrid,nslope), intent(inout) :: tendencies_h2o_ice_phys ! physical point field: Evolution of perennial ice over one year (kg/m^2/year) 40 logical, intent(inout) :: STOPPING ! Stopping criterion 41 41 42 42 ! local: 43 ! ---- 44 45 INTEGER :: i,j,islope ! loop variable 46 REAL :: pos_tend, neg_tend, real_coefficient,negative_part ! Variable to conserve water 47 REAL :: new_tendencies(ngrid,nslope) ! Tendencies computed in order to conserve water ice on the surface, only exchange between surface are done 48 43 ! ------ 44 integer :: i,j,islope ! loop variable 45 real :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve water 46 real :: new_tendencies(ngrid,nslope) ! Tendencies computed in order to conserve water ice on the surface, only exchange between surface are done 49 47 !======================================================================= 50 48 … … 53 51 pos_tend=0. 54 52 neg_tend=0. 55 if (ngrid.NE.1) then ! to make sure we are not in 1D 53 if (ngrid.NE.1) then ! to make sure we are not in 1D 56 54 ! We compute the amount of water accumulating and sublimating 57 55 do i=1,ngrid … … 77 75 enddo 78 76 ! We adapt the tendencies to conserve water and do only exchange between grid points 79 if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then ! We are sublimating more in the planet than condensing 77 if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then ! We are sublimating more in the planet than condensing 80 78 do i=1,ngrid 81 79 do islope=1,nslope … … 97 95 enddo 98 96 enddo 99 elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then 97 elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then 100 98 write(*,*) "Reason of stopping : There is either no water ice sublimating or no water ice increasing !!" 101 99 write(*,*) "Tendencies on ice sublimating=", neg_tend … … 123 121 enddo 124 122 enddo 125 126 123 124 127 125 if(pos_tend.eq.0) then 128 126 real_coefficient = 0. 129 else 130 real_coefficient = negative_part/pos_tend ! We compute a coefficient by which we should remove the ice that has been added 131 ! to places even if this ice was contributing to an unphysical negative amount 127 else 128 real_coefficient = negative_part/pos_tend ! We compute a coefficient by which we should remove the ice that has been added 129 ! to places even if this ice was contributing to an unphysical negative amount 132 130 ! of ice at other places 133 131 endif -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3129 r3130 42 42 use criterion_pem_stop_mod, only: criterion_waterice_stop, criterion_co2_stop 43 43 use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, & 44 m_noco2, threshold_water_frost2peren ial, threshold_co2_frost2perenial44 m_noco2, threshold_water_frost2perennial, threshold_co2_frost2perennial 45 45 use evol_co2_ice_s_mod, only: evol_co2_ice_s 46 46 use evol_h2o_ice_s_mod, only: evol_h2o_ice_s … … 77 77 albedodat, zmea, zstd, zsig, zgam, zthe, & 78 78 hmons, summit, base,albedo_h2o_frost, & 79 frost_albedo_threshold, emissiv, watercaptag, peren ial_co2ice, &79 frost_albedo_threshold, emissiv, watercaptag, perennial_co2ice, & 80 80 emisice, albedice 81 81 use dimradmars_mod, only: totcloudfrac, albedo … … 84 84 use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master 85 85 use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit 86 use paleoclimate_mod, only: albedo_peren ialco286 use paleoclimate_mod, only: albedo_perennialco2 87 87 use comcstfi_h, only: pi, rad, g, cpp, mugaz, r 88 88 #else … … 161 161 ! Variables for h2o_ice evolution 162 162 real :: ini_surf_h2o ! Initial surface of sublimating h2o ice 163 real :: ini_surf_co2 ! Initial surface of sublimating co2 ice164 163 real :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa] 165 164 real :: global_ave_press_old ! constant: Global average pressure of initial/previous time step … … 171 170 logical :: STOPPING_water ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached? 172 171 logical :: STOPPING_1_water ! Logical: is there still water ice to sublimate? 173 logical :: STOPPING_co2 ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached?174 172 logical :: STOPPING_pressure ! Logical: is the criterion (% of change in the surface pressure) reached? 175 173 integer :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param 176 174 real, save :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] 177 real, dimension(:,:), allocatable :: vmr_co2_gcm ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]178 real, dimension(:,:), allocatable :: vmr_co2_pem_phys ! Physics x Times co2 volume mixing ratio used in the PEM179 real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]180 175 real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg] 181 176 integer :: timelen ! # time samples … … 183 178 real, dimension(:,:), allocatable :: p ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) 184 179 real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 180 181 ! Variables for co2_ice evolution 182 real :: ini_surf_co2 ! Initial surface of sublimating co2 ice 183 logical :: STOPPING_co2 ! Logical: is the criterion (% of change in the surface of sublimating co2 ice) reached? 184 real, dimension(:,:), allocatable :: vmr_co2_gcm ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3] 185 real, dimension(:,:), allocatable :: vmr_co2_pem_phys ! Physics x Times co2 volume mixing ratio used in the PEM 186 real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg] 187 real, dimension(:,:), allocatable :: co2frost_ini ! Initial amount of co2 frost (at the start of the PEM run) 188 real, dimension(:,:), allocatable :: perennial_co2ice_ini ! Initial amoun of perennial co2 ice (at the start of the PEM run) 189 190 185 191 186 192 ! Variables for slopes … … 193 199 real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field: Logical array indicating if there is water ice at initial state 194 200 real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field: Logical array indicating if there is co2 ice at initial state 195 real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point x slope field: Tendency of evolution of peren ial co2 ice over a year196 real, dimension(:,:), allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of peren ial co2 ice over a year in the PCM197 real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical point x slope field: Tendency of evolution of peren ial h2o ice201 real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year 202 real, dimension(:,:), allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perennial co2 ice over a year in the PCM 203 real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical point x slope field: Tendency of evolution of perennial h2o ice 198 204 real, dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow 199 205 real, dimension(:), allocatable :: flag_co2flow_mesh ! (ngrid) : Flag where there is a CO2 glacier flow … … 383 389 call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, & 384 390 tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, & 385 watercap,peren ial_co2ice,def_slope,def_slope_mean,subslope_dist)391 watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist) 386 392 387 393 ! Remove unphysical values of surface tracer … … 588 594 if (tendencies_h2o_ice(i,islope) < 0) then 589 595 initial_h2o_ice(i,islope) = 1. 590 ini_surf_h2o =ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)596 ini_surf_h2o = ini_surf_h2o + cell_area(i)*subslope_dist(i,islope) 591 597 else 592 598 initial_h2o_ice(i,islope) = 0. … … 595 601 enddo 596 602 597 write(*,*) "Total initial surface of co2 ice sublimating (slope)=", ini_surf_co2598 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o603 write(*,*) "Total initial surface of co2 ice sublimating (slope) =", ini_surf_co2 604 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o 599 605 write(*,*) "Total surface of the planet=", Total_surface 600 606 allocate(zplev_gcm(ngrid,nlayer + 1)) … … 625 631 delta_h2o_icetablesublim(:) = 0. 626 632 633 ! We save the initial values for the co2 frost and perennial ice 634 allocate(co2frost_ini(ngrid,nslope),perennial_co2ice_ini(ngrid,nslope)) 635 co2frost_ini = qsurf(:,igcm_co2,:) 636 perennial_co2ice_ini = perennial_co2ice 637 627 638 do ig = 1,ngrid 628 639 do islope = 1,nslope 629 640 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.) 641 qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) + perennial_co2ice(ig,islope) 630 642 enddo 631 643 enddo … … 1000 1012 ! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap 1001 1013 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.)) 1002 qsurf(ig,igcm_h2o_ice,islope) =0.1014 qsurf(ig,igcm_h2o_ice,islope) = 0. 1003 1015 endif 1004 watercap_sum = watercap_sum +watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)1016 watercap_sum = watercap_sum + watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 1005 1017 watercap(ig,islope) = 0. 1006 1018 enddo … … 1015 1027 enddo 1016 1028 if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. 1017 if (water_sum > threshold_water_frost2peren ial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.11029 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 1018 1030 watercaptag(ig) = .true. 1019 water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2peren ial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost1031 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 1020 1032 do islope = 1,nslope 1021 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2peren ial/2.*cos(pi*def_slope_mean(islope)/180.)1033 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.) 1022 1034 enddo 1023 1035 endif 1024 1036 else ! let's check that the infinite source of water disapear 1025 if ((water_sum + water_reservoir(ig)) < threshold_water_frost2peren ial) then1037 if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perennial) then 1026 1038 watercaptag(ig) = .false. 1027 1039 do islope = 1,nslope … … 1036 1048 do ig = 1,ngrid 1037 1049 do islope = 1,nslope 1038 if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then 1039 perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope) 1040 qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope) 1041 albedo(ig,1,islope) = albedo_perenialco2 1042 albedo(ig,2,islope) = albedo_perenialco2 1050 if (qsurf(ig,igcm_co2,islope) == threshold_co2_frost2perennial) then 1051 ! If co2 ice is equal to the threshold, then everything is transformed into perennial ice 1052 perennial_co2ice(ig,islope) = threshold_co2_frost2perennial 1053 qsurf(ig,igcm_co2,islope) = 0. 1054 else if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perennial) then 1055 ! If co2 ice is superior to the threshold, then co2 frost is equal to the threshold (max possible value) 1056 ! and the leftover is transformed into perennial ice 1057 perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope) - threshold_co2_frost2perennial 1058 qsurf(ig,igcm_co2,islope) = threshold_co2_frost2perennial 1059 else 1060 ! If co2 ice is inferior to the threshold, then we compare with the initial state 1061 if (qsurf(ig,igcm_co2,islope) > perennial_co2ice_ini(ig,islope)) then 1062 ! If co2 ice higher than the initial perennial ice, then the change is affected only to the frost 1063 perennial_co2ice(ig,islope) = perennial_co2ice_ini(ig,islope) 1064 qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) - perennial_co2ice_ini(ig,islope) 1065 else 1066 ! If co2 ice is lower than the initial perennial ice, then there is no frost anymore 1067 perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope) 1068 qsurf(ig,igcm_co2,islope) = 0. 1069 endif 1043 1070 endif 1044 1071 enddo … … 1124 1151 nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, & 1125 1152 inertiedat,def_slope,subslope_dist) 1126 1127 1153 call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & 1128 1154 ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, & 1129 1155 albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, & 1130 wstar,watercap,peren ial_co2ice)1156 wstar,watercap,perennial_co2ice) 1131 1157 #else 1132 1158 call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, & 1133 1159 nlayer,nq,ptimestep,pday,time_phys,cell_area, & 1134 1160 albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) 1135 1136 1161 call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & 1137 1162 ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, & … … 1162 1187 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave) 1163 1188 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter) 1189 deallocate(co2frost_ini,perennial_co2ice_ini) 1164 1190 !----------------------------- END OUTPUT ------------------------------ 1165 1191 -
trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90
r3123 r3130 66 66 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen) :: h2o_ice_s_dyn ! h2o ice per slope of the concatenated file [kg/m^2] 67 67 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen) :: watercap 68 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen) :: perennial_co2ice 68 69 real, dimension(iim_input + 1,jjm_input + 1,timelen) :: vmr_co2_gcm ! CO2 volume mixing ratio in the first layer [mol/m^3] 69 70 real, dimension(iim_input + 1,jjm_input + 1,timelen) :: ps_GCM ! Surface Pressure [Pa] … … 106 107 call get_var3("watercap",watercap(:,:,1,:)) 107 108 write(*,*) "Data for watercap downloaded!" 109 110 call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:)) 111 write(*,*) "Data for perennial_co2ice downloaded!" 108 112 #endif 109 113 … … 142 146 enddo 143 147 write(*,*) "Data for watercap downloaded!" 148 149 do islope = 1,nslope 150 write(num,'(i2.2)') islope 151 call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:)) 152 enddo 153 write(*,*) "Data for perennial_co2ice downloaded!" 144 154 #endif 145 155 … … 179 189 write(*,*) "Computing the min of co2_ice..." 180 190 where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0. 181 min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn ,4)191 min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn + perennial_co2ice,4) 182 192 183 193 ! Compute averages over the year for each point … … 204 214 endif 205 215 if (q_h2o_dyn(i,j,t) < 0) then 206 q_h2o_dyn(i,j,t) = 1.e- 30216 q_h2o_dyn(i,j,t) = 1.e-10 207 217 else if (q_h2o_dyn(i,j,t) > 1) then 208 218 q_h2o_dyn(i,j,t) = 1. … … 282 292 !======================================================================= 283 293 284 SUBROUTINE get_var1( filename,v)285 286 implicit none 287 288 character(len = *), intent(in) :: filename294 SUBROUTINE get_var1(var,v) 295 296 implicit none 297 298 character(len = *), intent(in) :: var 289 299 real, dimension(:), intent(out) :: v 290 300 291 call error_msg(NF90_INQ_VARID(fID, filename,vID),"inq",filename)292 call error_msg(NF90_GET_VAR(fID,vID,v),"get", filename)301 call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var) 302 call error_msg(NF90_GET_VAR(fID,vID,v),"get",var) 293 303 294 304 END SUBROUTINE get_var1 … … 296 306 !======================================================================= 297 307 298 SUBROUTINE get_var3( filename,v) ! on U grid299 300 implicit none 301 302 character(len = *), intent(in) :: filename308 SUBROUTINE get_var3(var,v) ! on U grid 309 310 implicit none 311 312 character(len = *), intent(in) :: var 303 313 real, dimension(:,:,:), intent(out) :: v 304 314 305 call error_msg(NF90_INQ_VARID(fID, filename,vID),"inq",filename)306 call error_msg(NF90_GET_VAR(fID,vID,v),"get", filename)315 call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var) 316 call error_msg(NF90_GET_VAR(fID,vID,v),"get",var) 307 317 308 318 END SUBROUTINE get_var3 … … 310 320 !======================================================================= 311 321 312 SUBROUTINE get_var4( filename,v)313 314 implicit none 315 316 character(len = *), intent(in) :: filename322 SUBROUTINE get_var4(var,v) 323 324 implicit none 325 326 character(len = *), intent(in) :: var 317 327 real, dimension(:,:,:,:), intent(out) :: v 318 328 319 call error_msg(NF90_INQ_VARID(fID, filename,vID),"inq",filename)320 call error_msg(NF90_GET_VAR(fID,vID,v),"get", filename)329 call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var) 330 call error_msg(NF90_GET_VAR(fID,vID,v),"get",var) 321 331 322 332 END SUBROUTINE get_var4 -
trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90
r3076 r3130 29 29 real, intent(in) :: global_ave_press_GCM ! global averaged pressure at previous timestep 30 30 real, intent(in) :: global_ave_press_new ! global averaged pressure at current timestep 31 real, dimension(ngrid,nslope), intent(in) :: tendencies_co2_ice_phys_ini ! physical point field : Evolution of peren ial ice over one year31 real, dimension(ngrid,nslope), intent(in) :: tendencies_co2_ice_phys_ini ! physical point field : Evolution of perennial ice over one year 32 32 real, dimension(ngrid,nslope), intent(in) :: co2ice_slope ! CO2 ice per mesh and sub-grid slope(kg/m^2) 33 33 real, dimension(ngrid,nslope), intent(in) :: emissivity_slope ! Emissivity per mesh and sub-grid slope(1) 34 34 ! OUTPUT 35 real, intent(inout) :: tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of peren ial ice over one year35 real, intent(inout) :: tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perennial ice over one year 36 36 37 37 ! local: -
trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90
r2985 r3130 147 147 do islope=1,nslope 148 148 if (ice_depth(ig,islope).gt.-1e-6) then 149 ! 3.0 Case where it is peren ial ice149 ! 3.0 Case where it is perennial ice 150 150 if (ice_depth(ig,islope).lt. 1e-10) then 151 151 call ice_thermal_properties(.true.,1.,0.,ice_inertia)
Note: See TracChangeset
for help on using the changeset viewer.