Changeset 3493
- Timestamp:
- Nov 5, 2024, 5:19:11 PM (4 weeks ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 1 deleted
- 5 edited
- 15 moved
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90
r3492 r3493 12 12 13 13 14 SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,T _in,nsoil,thIn,p0,pfrost,porefill_in,porefill,ssi_depth)14 SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,Tb,nz,thIn,p0,pfrost,porefill_in,porefill,ssi_depth) 15 15 16 16 !*********************************************************************** … … 22 22 implicit none 23 23 integer, parameter :: NP=1 ! # of sites 24 integer nz, i, k, nsoil,iloop24 integer nz, i, k, iloop 25 25 real(8) zmax, delta, z(NMAX), icetime, porosity, icefrac 26 26 real(8), dimension(NP) :: albedo, thIn, rhoc … … 30 30 real(8), dimension(NP) :: zdepthF, zdepthE, zdepthT, zdepthG 31 31 real(8), dimension(NMAX,NP) :: porefill, porefill_in 32 real(8), dimension(NP) :: Tb, T_in, Tmean1, Tmean3, avrho1 32 real(8), dimension(nz) :: Tb 33 real(8), dimension(NP) :: Tmean1, Tmean3, avrho1 33 34 real(8) tmax, tlast, avrho1prescribed(NP), l1 34 35 real(8), external :: smartzfac … … 42 43 43 44 ! parameters that never ever change 44 nz=nsoil45 45 porosity = 0.4d0 ! porosity of till 46 46 !rhoc(:) = 1500.*800. ! will be overwritten … … 69 69 !call setgrid(nz,z,zmax,zfac) 70 70 l1=2.e-4 71 do iloop=0,n soil71 do iloop=0,nz 72 72 z(iloop) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.))) 73 73 enddo … … 89 89 90 90 ! initializations 91 Tb=T_in91 !Tb = -9999. 92 92 zdepthF(:) = -9999. 93 93 -
trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90
r3492 r3493 32 32 real(8), intent(IN) :: bigstep 33 33 real(8), intent(IN) :: thIn(NP), rhoc(NP), z(NMAX), porosity, pfrost(NP) 34 real(8), intent(INOUT) :: Tb(NP), porefill(nz,NP), zdepthF(NP), zdepthT(NP) 34 real(8), intent(INOUT) :: porefill(nz,NP), zdepthF(NP), zdepthT(NP) 35 real(8), intent(INOUT) :: Tb(nz) 35 36 real(8), intent(OUT), dimension(NP) :: zdepthE, Tmean1, Tmean3, zdepthG 36 37 real(8), intent(IN), dimension(NP) :: albedo, p0 … … 149 150 real(8), intent(OUT) :: avdrho, avdrhoP ! difference in annual mean vapor density 150 151 real(8), intent(OUT) :: avrho1 ! mean vapor density on surface 151 real(8), intent(INOUT) :: Tb, Tmean1 152 real(8), intent(INOUT) :: Tmean1 153 real(8), intent(INOUT) :: Tb(nz) 152 154 integer, intent(OUT) :: typeF ! index of depth below which filling occurs 153 155 real(8), intent(OUT) :: zdepthE, zdepthF -
trunk/LMDZ.COMMON/libf/evolution/NS_generalorbit.F90
r3492 r3493 1 C=============================================================2 CSubroutine to return distance, longitude, and declination of3 Cthe sun in planetocentric coordinates from orbital elements4 C 5 CINPUTS:6 Cedays = time since perihelion (earth days)7 Ca = semimajor axis (AU)8 Cecc = eccentricity9 Comega = Ls of perihelion, relative to equinox (radians)10 Ceps = obliquity (radians)11 C 12 COUTPUTS:13 CLs = areocentric longitude (radians)14 Cdec = planetocentric solar declination (radians)15 Cr = heliocentric distance (AU)16 C=============================================================1 !============================================================= 2 ! Subroutine to return distance, longitude, and declination of 3 ! the sun in planetocentric coordinates from orbital elements 4 ! 5 ! INPUTS: 6 ! edays = time since perihelion (earth days) 7 ! a = semimajor axis (AU) 8 ! ecc = eccentricity 9 ! omega = Ls of perihelion, relative to equinox (radians) 10 ! eps = obliquity (radians) 11 ! 12 ! OUTPUTS: 13 ! Ls = areocentric longitude (radians) 14 ! dec = planetocentric solar declination (radians) 15 ! r = heliocentric distance (AU) 16 !============================================================= 17 17 18 18 SUBROUTINE generalorbit(edays,a,ecc,omega,eps,Ls,dec,r) … … 25 25 real*8 M,E,nu,T,Eold 26 26 27 cT = orbital period (days)27 ! T = orbital period (days) 28 28 T = sqrt(4*pi**2/(6.674e-11*1.989e30)*(a*149.598e9)**3)/86400. 29 29 30 cM = mean anomaly (radians)30 ! M = mean anomaly (radians) 31 31 M = 2.*pi*edays/T ! M=0 at perihelion 32 32 33 cE = eccentric anomaly34 csolve M = E - ecc*sin(E) by newton method33 ! E = eccentric anomaly 34 ! solve M = E - ecc*sin(E) by newton method 35 35 E = M 36 36 do j=1,10 … … 40 40 enddo 41 41 42 cnu = true anomaly42 ! nu = true anomaly 43 43 !nu = acos(cos(E)-ecc/(1.-ecc*cos(E))) 44 44 !nu = sqrt(1-ecc^2)*sin(E)/(1.-ecc*cos(E)) -
trunk/LMDZ.COMMON/libf/evolution/NS_psv.F90
r3492 r3493 1 1 real*8 function psv(T) 2 Csaturation vapor pressure of H2O ice [Pascal]3 Cinput is temperature [Kelvin]2 ! saturation vapor pressure of H2O ice [Pascal] 3 ! input is temperature [Kelvin] 4 4 implicit none 5 5 real*8 T 6 6 7 C-----parametrization 18 creal*8 DHmelt,DHvap,DHsub,R,pt,Tt,C9 cparameter (DHmelt=6008.,DHvap=45050.)10 cparameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]11 cparameter (R=8.314,pt=6.11e2,Tt=273.16)12 cC = (DHsub/R)*(1./T - 1./Tt)13 cpsv = pt*exp(-C)7 !-----parametrization 1 8 ! real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C 9 ! parameter (DHmelt=6008.,DHvap=45050.) 10 ! parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol] 11 ! parameter (R=8.314,pt=6.11e2,Tt=273.16) 12 ! C = (DHsub/R)*(1./T - 1./Tt) 13 ! psv = pt*exp(-C) 14 14 15 C-----parametrization 216 Ceq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)17 Cdiffers from parametrization 1 by only 0.1%15 !-----parametrization 2 16 ! eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005) 17 ! differs from parametrization 1 by only 0.1% 18 18 real*8 A,B 19 19 parameter (A=-6143.7, B=28.9074) 20 20 psv = exp(A/T+B) ! Clapeyron 21 21 22 C-----parametrization 323 Ceq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)24 cpsv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)22 !-----parametrization 3 23 ! eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005) 24 ! psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T) 25 25 26 26 end … … 29 29 30 30 real*8 function frostpoint(p) 31 Cinverse of psv32 Cinput is partial pressure [Pascal]33 Coutput is temperature [Kelvin]31 ! inverse of psv 32 ! input is partial pressure [Pascal] 33 ! output is temperature [Kelvin] 34 34 implicit none 35 35 real*8 p 36 36 37 C-----inverse of parametrization 138 creal*8 DHmelt,DHvap,DHsub,R,pt,Tt39 cparameter (DHmelt=6008.,DHvap=45050.)40 cparameter (DHsub=DHmelt+DHvap)41 cparameter (R=8.314,pt=6.11e2,Tt=273.16)42 cfrostpoint = 1./(1./Tt-R/DHsub*log(p/pt))37 !-----inverse of parametrization 1 38 ! real*8 DHmelt,DHvap,DHsub,R,pt,Tt 39 ! parameter (DHmelt=6008.,DHvap=45050.) 40 ! parameter (DHsub=DHmelt+DHvap) 41 ! parameter (R=8.314,pt=6.11e2,Tt=273.16) 42 ! frostpoint = 1./(1./Tt-R/DHsub*log(p/pt)) 43 43 44 C-----inverse of parametrization 245 Cinverse of eq. (2) in Murphy & Koop (2005)44 !-----inverse of parametrization 2 45 ! inverse of eq. (2) in Murphy & Koop (2005) 46 46 real*8 A,B 47 47 parameter (A=-6143.7, B=28.9074) 48 48 frostpoint = A / (log(p) - B) 49 49 50 C-----approximate inverse of parametrization 351 Ceq. (8) in Murphy & Koop (2005)52 cfrostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))50 !-----approximate inverse of parametrization 3 51 ! eq. (8) in Murphy & Koop (2005) 52 ! frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p)) 53 53 54 54 end -
trunk/LMDZ.COMMON/libf/evolution/NS_psvco2.F90
r3492 r3493 1 1 real*8 function psvco2(T) 2 Csolid-vapor transition for CO23 Creturns saturation pressure in Pascal, input is temperature in Kelvin2 ! solid-vapor transition for CO2 3 ! returns saturation pressure in Pascal, input is temperature in Kelvin 4 4 implicit none 5 5 real*8 T … … 10 10 11 11 real*8 function tfrostco2(p) 12 Cthe inverse of function psvco213 Cinput is pressure in Pascal, output is temperature in Kelvin12 ! the inverse of function psvco2 13 ! input is pressure in Pascal, output is temperature in Kelvin 14 14 implicit none 15 15 real*8 p … … 20 20 21 21 22 C------------------------------------------------------------------------------23 CAntoine equation parameters from NIST, 154K-196K24 Cbased on Giauque and Egan (1937)25 CA=6.81228, B=1301.679, C=-3.49426 cp = 10**(A-(B/(T+C)))*1.e522 !------------------------------------------------------------------------------ 23 ! Antoine equation parameters from NIST, 154K-196K 24 ! based on Giauque and Egan (1937) 25 ! A=6.81228, B=1301.679, C=-3.494 26 ! p = 10**(A-(B/(T+C)))*1.e5 27 27 28 CExpressions from Int. Crit. Tabl. 3.207, based on many references29 Cmm2Pa = 133.3230 C-135C to -56.7C (138K-216K)31 cp = 10**(-0.05223*26179.3/T + 9.9082)*mm2Pa32 C-183C to -135C (90K-138K)33 cp = 10**(-1275.62/T + 0.006833*T + 8.3071)*mm2Pa34 CExpressions from Int. Crit. Tabl. 3.208, based on Henning35 C-110 to -80C (163K-193K)36 cp = 10**(- 1279.11/T + 1.75*log10(T) - 0.0020757*T + 5.85242)*mm2Pa37 cp = 10**(- 1352.6/T + 9.8318)*mm2Pa28 ! Expressions from Int. Crit. Tabl. 3.207, based on many references 29 ! mm2Pa = 133.32 30 ! -135C to -56.7C (138K-216K) 31 ! p = 10**(-0.05223*26179.3/T + 9.9082)*mm2Pa 32 ! -183C to -135C (90K-138K) 33 ! p = 10**(-1275.62/T + 0.006833*T + 8.3071)*mm2Pa 34 ! Expressions from Int. Crit. Tabl. 3.208, based on Henning 35 ! -110 to -80C (163K-193K) 36 ! p = 10**(- 1279.11/T + 1.75*log10(T) - 0.0020757*T + 5.85242)*mm2Pa 37 ! p = 10**(- 1352.6/T + 9.8318)*mm2Pa 38 38 39 CMars book (1992), p959, fit by chapter authors40 cp = exp(23.3494 - 3182.48/T)*100. ! 120-160K41 cp = exp(25.2194 - 3311.57/T - 6.71e-3*T)*100 ! 100-170K42 CMars book (1992), p960, based on Miller & Smythe (1970)43 cp = exp(26.1228 - 3385.26/T - 9.4461e-3*T)*100 ! 123-173K39 ! Mars book (1992), p959, fit by chapter authors 40 ! p = exp(23.3494 - 3182.48/T)*100. ! 120-160K 41 ! p = exp(25.2194 - 3311.57/T - 6.71e-3*T)*100 ! 100-170K 42 ! Mars book (1992), p960, based on Miller & Smythe (1970) 43 ! p = exp(26.1228 - 3385.26/T - 9.4461e-3*T)*100 ! 123-173K 44 44 45 CFray & Schmitt, PSS 57, 2053 (2009)46 CA0=1.476e1, A1=-2.571e3, A2=-7.781e4, A3=4.325e6, A4=-1.207e8, A5=1.350e947 cp = exp(A0+A1/T+A2/T**2+A3/T**3+A4/T**4+A5/T**5)*1e5 ! 40K-194.7K48 cA0=1.861e1, A1=-4.154e3, A2=1.041e549 cp = exp(A0 + A1/T + A2/T**2)*1e5 ! 194.7K-216.58K50 C------------------------------------------------------------------------------45 ! Fray & Schmitt, PSS 57, 2053 (2009) 46 ! A0=1.476e1, A1=-2.571e3, A2=-7.781e4, A3=4.325e6, A4=-1.207e8, A5=1.350e9 47 ! p = exp(A0+A1/T+A2/T**2+A3/T**3+A4/T**4+A5/T**5)*1e5 ! 40K-194.7K 48 ! A0=1.861e1, A1=-4.154e3, A2=1.041e5 49 ! p = exp(A0 + A1/T + A2/T**2)*1e5 ! 194.7K-216.58K 50 !------------------------------------------------------------------------------ -
trunk/LMDZ.COMMON/libf/evolution/NS_tridag.F90
r3492 r3493 1 C================================================2 CTridiagonal solver3 C================================================1 !================================================ 2 ! Tridiagonal solver 3 !================================================ 4 4 SUBROUTINE tridag(a,b,c,r,u,n) 5 5 INTEGER n,NMAX … … 11 11 stop 'tridag: rewrite equations' 12 12 endif 13 cif(n.gt.NMAX) then14 cprint *, 'tridag: too many points, set NMAX>',n15 cstop16 cendif13 ! if(n.gt.NMAX) then 14 ! print *, 'tridag: too many points, set NMAX>',n 15 ! stop 16 ! endif 17 17 bet=b(1) 18 18 u(1)=r(1)/bet 19 do 11j=2,n19 do j=2,n 20 20 gam(j)=c(j-1)/bet 21 21 bet=b(j)-a(j)*gam(j) 22 cif(bet.eq.0.)pause 'tridag failed'22 ! if(bet.eq.0.)pause 'tridag failed' 23 23 u(j)=(r(j)-a(j)*u(j-1))/bet 24 11 continue 25 do 12j=n-1,1,-124 enddo 25 do j=n-1,1,-1 26 26 u(j)=u(j)-gam(j+1)*u(j+1) 27 12 continue 28 return27 enddo 28 ! return 29 29 END 30 C(C) Copr. 1986-92 Numerical Recipes Software 0(9p#31(+.30 ! (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31(+. -
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3492 r3493 460 460 == 05/11/2024 == JBC 461 461 Correction in the formula to update the potential temperature introduced in r3442. 462 463 == 05/11/2024 == JBC 464 - Renaming of Norbert Schorghofer's subroutines with the prefix 'NS_'; 465 - Making the extension of all NS's subroutines as '.F90'; 466 - Deletion of the wrapper subroutine; 467 - Making the initialization, variables management and arguments of the main subroutine for the dynamic computation of ice table to be more suitable. -
trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
r3490 r3493 10 10 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 11 11 12 logical, save :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium 13 logical, save :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method 14 real, allocatable, dimension(:,:), save :: porefillingice_depth ! ngrid x nslope: Depth of the ice table [m] 15 real, allocatable, dimension(:,:), save :: porefillingice_thickness ! ngrid x nslope: Thickness of the ice table [m] 12 logical, save :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium 13 logical, save :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method 14 real, allocatable, dimension(:,:) :: icetable_depth ! ngrid x nslope: Depth of the ice table [m] 15 real, allocatable, dimension(:,:) :: icetable_thickness ! ngrid x nslope: Thickness of the ice table [m] 16 real, allocatable, dimension(:,:,:) :: ice_porefilling ! the amout of porefilling in each layer in each grid [m^3/m^3] 16 17 17 18 !----------------------------------------------------------------------- 18 19 contains 19 20 !----------------------------------------------------------------------- 20 SUBROUTINE ini_ice_table _porefilling(ngrid,nslope)21 SUBROUTINE ini_ice_table(ngrid,nslope,nsoil) 21 22 22 23 implicit none … … 24 25 integer, intent(in) :: ngrid ! number of atmospheric columns 25 26 integer, intent(in) :: nslope ! number of slope within a mesh 26 27 allocate(porefillingice_depth(ngrid,nslope)) 28 allocate(porefillingice_thickness(ngrid,nslope)) 29 30 END SUBROUTINE ini_ice_table_porefilling 31 32 !----------------------------------------------------------------------- 33 SUBROUTINE end_ice_table_porefilling 34 35 implicit none 36 37 if (allocated(porefillingice_depth)) deallocate(porefillingice_depth) 38 if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness) 39 40 END SUBROUTINE end_ice_table_porefilling 27 integer, intent(in) :: nsoil ! number of soil layers 28 29 allocate(icetable_depth(ngrid,nslope)) 30 if (icetable_equilibrium) then 31 allocate(icetable_thickness(ngrid,nslope)) 32 else if (icetable_dynamic) then 33 allocate(ice_porefilling(ngrid,nsoil,nslope)) 34 endif 35 36 END SUBROUTINE ini_ice_table 37 38 !----------------------------------------------------------------------- 39 SUBROUTINE end_ice_table 40 41 implicit none 42 43 if (allocated(icetable_depth)) deallocate(icetable_depth) 44 if (allocated(icetable_thickness)) deallocate(icetable_thickness) 45 if (allocated(ice_porefilling)) deallocate(ice_porefilling) 46 47 END SUBROUTINE end_ice_table 41 48 42 49 !------------------------------------------------------------------------------------------------------ -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3492 r3493 54 54 use orbit_param_criterion_mod, only: orbit_param_criterion 55 55 use recomp_orb_param_mod, only: recomp_orb_param 56 use ice_table_mod, only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &57 ini_ice_table _porefilling, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium,compute_massh2o_exchange_ssi56 use ice_table_mod, only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, & 57 ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi 58 58 use soil_thermalproperties_mod, only: update_soil_thermalproperties 59 59 use time_phylmdz_mod, only: daysec, dtphys … … 204 204 205 205 ! Variables for surface and soil 206 real, dimension(:,:), allocatable :: tsurf_avg ! Physic x SLOPE field: Averaged Surface Temperature [K] 207 real, dimension(:,:,:), allocatable :: tsoil_avg ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K] 208 real, dimension(:,:,:), allocatable :: tsurf_PCM_timeseries ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K] 209 real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K] 210 real, dimension(:,:,:,:), allocatable :: tsoil_anom ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K] 211 real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K] 212 real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: Intermediate when computing Tsoil [K] 213 real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] 214 real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] 215 real, dimension(:,:), allocatable :: watersurf_density_avg ! Physic x Slope, water surface density, yearly averaged [kg/m^3] 216 real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] 217 real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3] 218 real, dimension(:,:), allocatable :: Tsurfavg_before_saved ! Surface temperature saved from previous time step [K] 219 real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] 220 real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] 221 real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] 222 real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] 223 logical :: bool_sublim ! logical to check if there is sublimation or not 224 logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep 225 real, dimension(:,:), allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] 226 real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] 227 228 ! Dynamic ice table 229 real, dimension(:,:), allocatable :: ss_ice_pf ! the amout of porefilling in each layer in each grid [m^3/m^3] 230 real, dimension(:,:), allocatable :: ss_ice_pf_new ! the amout of porefilling in each layer in each grid after the ss module[m^3/m^3] 231 real, dimension(:,:), allocatable :: porefillingice_depth_new ! new pure ice table depth 206 real, dimension(:,:), allocatable :: tsurf_avg ! Physic x SLOPE field: Averaged Surface Temperature [K] 207 real, dimension(:,:,:), allocatable :: tsoil_avg ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K] 208 real, dimension(:,:,:), allocatable :: tsurf_PCM_timeseries ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K] 209 real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K] 210 real, dimension(:,:,:,:), allocatable :: tsoil_anom ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K] 211 real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K] 212 real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: Intermediate when computing Tsoil [K] 213 real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] 214 real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] 215 real, dimension(:,:), allocatable :: watersurf_density_avg ! Physic x Slope, water surface density, yearly averaged [kg/m^3] 216 real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] 217 real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3] 218 real, dimension(:,:), allocatable :: Tsurfavg_before_saved ! Surface temperature saved from previous time step [K] 219 real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] 220 real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] 221 real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] 222 real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] 223 logical :: bool_sublim ! logical to check if there is sublimation or not 224 logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep 225 real, dimension(:,:), allocatable :: icetable_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] 226 real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] 227 real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table 228 real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table 229 232 230 ! Some variables for the PEM run 233 231 real, parameter :: year_step = 1 ! Timestep for the pem … … 557 555 allocate(delta_co2_adsorbded(ngrid)) 558 556 allocate(co2ice_disappeared(ngrid,nslope)) 559 allocate( porefillingice_thickness_prev_iter(ngrid,nslope))557 allocate(icetable_thickness_prev_iter(ngrid,nslope)) 560 558 allocate(delta_h2o_icetablesublim(ngrid)) 561 559 allocate(delta_h2o_adsorbded(ngrid)) … … 583 581 call end_adsorption_h_PEM 584 582 call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) 585 call end_ice_table _porefilling586 call ini_ice_table _porefilling(ngrid,nslope)583 call end_ice_table 584 call ini_ice_table(ngrid,nslope,nsoilmx_PEM) 587 585 588 586 if (soil_pem) then … … 646 644 endif 647 645 648 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM, porefillingice_depth, &649 porefillingice_thickness,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,&650 tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,global_avg_press_PCM,&651 watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,delta_co2_adsorbded,&652 h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)646 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, & 647 icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys, & 648 ps_timeseries,tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, & 649 global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys, & 650 delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif) 653 651 654 652 ! We save the places where h2o ice is sublimating … … 928 926 ! II_d.3 Update the ice table 929 927 if (icetable_equilibrium) then 930 write(*,*) "Compute ice table "931 porefillingice_thickness_prev_iter = porefillingice_thickness932 call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:), porefillingice_depth,porefillingice_thickness)933 call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM, porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere928 write(*,*) "Compute ice table (equilibrium method)" 929 icetable_thickness_prev_iter = icetable_thickness 930 call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness) 931 call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere 934 932 else if (icetable_dynamic) then 935 write(*,*) "Compute ice table" 936 porefillingice_thickness_prev_iter = porefillingice_thickness 937 call dyn_ss_ice_m_wrapper(ngrid,nsoilmx,TI_PEM,ps,mmol(igcm_h2o_vap),tsoil_PEM,porefillingice_depth,ss_ice_pf,ss_ice_pf_new,porefillingice_depth_new) 938 939 call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere 933 write(*,*) "Compute ice table (dynamic method)" 934 allocate(porefill(nsoilmx_PEM)) 935 do ig = 1,ngrid 936 do islope = 1,nslope 937 call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2),ice_porefilling(ig,:,islope),porefill,ssi_depth) 938 icetable_depth(ig,islope) = ssi_depth 939 ice_porefilling(ig,:,islope) = porefill 940 enddo 941 enddo 942 deallocate(porefill) 943 call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere 940 944 endif 941 945 942 946 ! II_d.4 Update the soil thermal properties 943 call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new, porefillingice_depth,porefillingice_thickness,TI_PEM)947 call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM) 944 948 945 949 ! II_d.5 Update the mass of the regolith adsorbed … … 952 956 totmassh2o_adsorbded = 0. 953 957 do ig = 1,ngrid 954 do islope = 1,nslope958 do islope = 1,nslope 955 959 do l = 1,nsoilmx_PEM 956 if (l ==1) then960 if (l == 1) then 957 961 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & 958 962 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) … … 960 964 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 961 965 else 962 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l -1))* &966 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & 963 967 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 964 totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l -1))* &968 totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & 965 969 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 966 970 endif … … 987 991 call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) 988 992 if (icetable_equilibrium) then 989 call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2, porefillingice_depth(:,islope))990 call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2, porefillingice_thickness(:,islope))993 call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) 994 call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope)) 991 995 else if (icetable_dynamic) then 992 call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2, porefillingice_depth(:,islope))993 call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2, porefillingice_thickness(:,islope))996 call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) 997 call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope)) 994 998 endif 995 999 … … 1205 1209 call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist) 1206 1210 call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, & 1207 TI_PEM, porefillingice_depth,porefillingice_thickness,&1211 TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, & 1208 1212 co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif) 1209 1213 write(*,*) "restartpem.nc has been written" … … 1227 1231 deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved) 1228 1232 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg) 1229 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim, porefillingice_thickness_prev_iter)1233 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,icetable_thickness_prev_iter) 1230 1234 deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif) 1231 1235 !----------------------------- END OUTPUT ------------------------------ -
trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
r3490 r3493 7 7 !======================================================================= 8 8 9 SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table ,ice_table_thickness,&10 tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,&11 global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys, &9 SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness, & 10 ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, & 11 global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys, & 12 12 m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif) 13 13 … … 60 60 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! soil (mid-layer) thermal inertia in the PEM grid [SI] 61 61 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM ! soil (mid-layer) temperature [K] 62 real, dimension(ngrid,nslope), intent(inout) :: ice_table 62 real, dimension(ngrid,nslope), intent(inout) :: ice_table_depth ! Ice table depth [m] 63 63 real, dimension(ngrid,nslope), intent(inout) :: ice_table_thickness ! Ice table thickness [m] 64 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling ! Subsurface ice pore filling [m3/m3] 64 65 real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst ! instantaneous soil (mid-layer) temperature [K] 65 66 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2] … … 107 108 108 109 !0.2 Set to default values 109 ice_table = -1.! by default, no ice table110 ice_table_depth = -1. ! by default, no ice table 110 111 ice_table_thickness = -1. 111 112 !1. Run … … 225 226 do ig = 1,ngrid 226 227 if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then 227 inertiedat_PEM(ig,index_breccia +1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &228 (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &228 inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ & 229 (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + & 229 230 ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2)))) 230 231 … … 308 309 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 309 310 !3. Ice Table 310 if (icetable_equilibrium .or. icetable_dynamic) then311 call get_field("ice_table ",ice_table,found)311 if (icetable_equilibrium) then 312 call get_field("ice_table_depth",ice_table_depth,found) 312 313 if (.not. found) then 313 write(*,*)'PEM settings: failed loading <ice_table >'314 write(*,*)'PEM settings: failed loading <ice_table_depth>' 314 315 write(*,*)'will reconstruct the values of the ice table given the current state' 315 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table ,ice_table_thickness)316 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table ,ice_table_thickness,TI_PEM)316 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness) 317 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM) 317 318 do islope = 1,nslope 318 319 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) 319 320 enddo 321 endif 322 write(*,*) 'PEMETAT0: ICE TABLE done' 323 else if (icetable_dynamic) then 324 call get_field("ice_table_depth",ice_table_depth,found) 325 if (.not. found) then 326 write(*,*)'PEM settings: failed loading <ice_table_depth>' 327 write(*,*)'will reconstruct the values of the ice table given the current state' 328 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness) 329 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM) 330 do islope = 1,nslope 331 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) 332 enddo 333 endif 334 call get_field("ice_porefilling",ice_porefilling,found) 335 if (.not. found) then 336 write(*,*)'PEM settings: failed loading <ice_porefilling>' 337 ice_porefilling = 0. 320 338 endif 321 339 write(*,*) 'PEMETAT0: ICE TABLE done' … … 354 372 call close_startphy 355 373 356 else !No start fi, let's build all by hand374 else !No startpem, let's build all by hand 357 375 358 376 write(*,*)'There is no "startpem.nc"!' … … 484 502 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 485 503 !c) Ice table 486 if (icetable_equilibrium .or. icetable_dynamic) then 487 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness) 488 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM) 504 if (icetable_equilibrium) then 505 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness) 506 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM) 507 do islope = 1,nslope 508 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) 509 enddo 510 write(*,*) 'PEMETAT0: Ice table done' 511 else if (icetable_dynamic) then 512 ice_porefilling = 0. 513 ice_table_depth = 0. 514 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM) 489 515 do islope = 1,nslope 490 516 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) -
trunk/LMDZ.COMMON/libf/evolution/pemredem.F90
r3319 r3493 58 58 59 59 SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, & 60 ice_table_depth,ice_table_thickness, m_co2_regolith,m_h2o_regolith,h2o_ice,stratif)60 ice_table_depth,ice_table_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif) 61 61 62 62 ! write time-dependent variable to restart file … … 78 78 real, dimension(ngrid,nslope), intent(in) :: ice_table_depth ! under mesh bining according to slope 79 79 real, dimension(ngrid,nslope), intent(in) :: ice_table_thickness ! under mesh bining according to slope 80 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: ice_porefilling ! under mesh bining according to slope 80 81 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith 81 82 real, dimension(ngrid,nslope), intent(in) :: h2o_ice … … 122 123 call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year) 123 124 call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year) 125 call put_field("ice_porefilling","Subsurface ice pore filling",ice_porefilling,Year) 124 126 call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year) 125 127 endif ! soil_pem
Note: See TracChangeset
for help on using the changeset viewer.