Index: trunk/LMDZ.COMMON/libf/evolution/abort_pem.F =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/abort_pem.F (revision 2993) +++ trunk/LMDZ.COMMON/libf/evolution/abort_pem.F (revision 2995) @@ -32,5 +32,5 @@ character(len=*), intent(in):: message - write(lunout,*) 'in abort_gcm' + write(lunout,*) 'in abort_pem' #ifdef CPP_XIOS Index: trunk/LMDZ.COMMON/libf/evolution/co2glaciers_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/co2glaciers_mod.F90 (revision 2993) +++ (revision ) @@ -1,240 +1,0 @@ -module co2glaciers_mod - - implicit none - LOGICAL co2glaciersflow ! True by default, to compute co2 ice flow. Read in pem.def - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! -!!! Purpose: Compute CO2 glacier flows -!!! -!!! Author: LL -!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -contains - - -subroutine co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,co2ice_slope,flag_co2flow,flag_co2flow_mesh) - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! -!!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do -!!! the ice transfer -!!! -!!! -!!! Author: LL -!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - -IMPLICIT NONE - -! arguments -! --------- - -! Inputs: - INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! number of time sample, physical points, subslopes, index of the flat subslope - REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles - REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol] - REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical x Time field: surface pressure given by the GCM [Pa] - REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surface pressure from the GCM [Pa] - REAL,INTENT(IN) :: global_ave_ps_PEM ! global averaged surface pressure during the PEM iteration [Pa] - -! Ouputs: - REAL,INTENT(INOUT) :: co2ice_slope(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] - REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes - REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid) ! same but within the mesh - - -! Local - REAL :: Tcond(ngrid) ! Physical field: CO2 condensation temperature [K] - REAL :: hmax(ngrid,nslope) ! Physical x Slope field: maximum thickness for co2 glacier before flow - -!----------------------------- - call computeTcond(timelen,ngrid,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond) - - call compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax) - - call transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh) - RETURN -end subroutine - - - -subroutine compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax) - - USE comconst_mod, ONLY: pi,g - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! -!!! Purpose: Compute the maximum thickness of CO2 glaciers given a slope angle -!!! before initating flow -!!! -!!! Author: LL,based on work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022) -!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - IMPLICIT NONE - -! arguments -! -------- - -! Inputs - INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes - INTEGER,INTENT(IN) :: iflat ! index of the flat subslope - REAL,INTENT(IN) :: Tcond(ngrid) ! Physical field: CO2 condensation temperature [K] - REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg] -! Outputs - REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum co2 thickness before flaw [m] -! Local - DOUBLE PRECISION :: tau_d = 5.e3 ! characteristic basal drag, understood as the stress that an ice CO2 mass flowing under its weight balanced by viscosity. Value obtained from I.Smith - REAL :: rho_co2(ngrid) ! co2 ice density [kg/m^3] - INTEGER :: ig,islope ! loop variables - REAL :: slo_angle - -! 1. Compute rho_co2 - DO ig = 1,ngrid - rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3 ! Mangan et al. 2017 - ENDDO - -! 3. Compute max thickness - DO ig = 1,ngrid - DO islope = 1,nslope - if(islope.eq.iflat) then - hmax(ig,islope) = 1.e6 - else - slo_angle = abs(def_slope_mean(islope)*pi/180.) - hmax(ig,islope) = tau_d/(rho_co2(ig)*g*slo_angle) - endif - ENDDO - ENDDO -RETURN - -end subroutine - -subroutine transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! -!!! Purpose: Transfer the excess of ice from one subslope to another -!!! No transfer between mesh at the time -!!! Author: LL -!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - USE comconst_mod, ONLY: pi - - -implicit none - -! arguments -! -------- - -! Inputs - INTEGER, INTENT(IN) :: ngrid,nslope !# of physical points and subslope - INTEGER, INTENT(IN) :: iflat ! index of the flat subslope - REAL, INTENT(IN) :: subslope_dist(ngrid,nslope) ! Distribution of the subgrid slopes within the mesh - REAL, INTENT(IN) :: def_slope_mean(nslope) ! values of the subgrid slopes - REAL, INTENT(IN) :: hmax(ngrid,nslope) ! maximum height of the CO2 glaciers before initiating flow [m] - REAL, INTENT(IN) :: Tcond(ngrid) ! CO2 condensation temperature [K] -! Outputs - REAL, INTENT(INOUT) :: co2ice_slope(ngrid,nslope) ! CO2 in the subslope [kg/m^2] - REAL, INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope - REAL, INTENT(INOUT) :: flag_co2flow_mesh(ngrid) ! boolean to check if there is flow in the mesh -! Local - INTEGER ig,islope ! loop - REAL rho_co2(ngrid) ! density of CO2, temperature dependant [kg/m^3] - INTEGER iaval ! ice will be transfered here - - - -! 0. Compute rho_co2 - DO ig = 1,ngrid - rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3 ! Mangan et al. 2017 - ENDDO - -! 1. Compute the transfer of ice - - DO ig = 1,ngrid - DO islope = 1,nslope - IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground -! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely on flat ground - IF(co2ice_slope(ig,islope).ge.rho_co2(ig)*hmax(ig,islope) * & - cos(pi*def_slope_mean(islope)/180.)) THEN -! Second: determine the flatest slopes possible: - IF(islope.gt.iflat) THEN - iaval=islope-1 - ELSE - iaval=islope+1 - ENDIF - do while ((iaval.ne.iflat).and. & - (subslope_dist(ig,iaval).eq.0)) - IF(iaval.gt.iflat) THEN - iaval=iaval-1 - ELSE - iaval=iaval+1 - ENDIF - enddo - co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + & - (co2ice_slope(ig,islope) - rho_co2(ig)*hmax(ig,islope) * & - cos(pi*def_slope_mean(islope)/180.)) * & - subslope_dist(ig,islope)/subslope_dist(ig,iaval) * & - cos(pi*def_slope_mean(iaval)/180.) / & - cos(pi*def_slope_mean(islope)/180.) - - co2ice_slope(ig,islope)=rho_co2(ig)*hmax(ig,islope) * & - cos(pi*def_slope_mean(islope)/180.) - - flag_co2flow(ig,islope) = 1. - flag_co2flow_mesh(ig) = 1. - ENDIF ! co2ice > hmax - ENDIF ! iflat - ENDDO !islope - ENDDO !ig -RETURN -end subroutine - - subroutine computeTcond(timelen,ngrid,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! -!!! Purpose: Compute CO2 condensation temperature -!!! -!!! Author: LL -!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2 - -implicit none - -! arguments: -! ---------- - -! INPUT - INTEGER,INTENT(IN) :: timelen, ngrid ! # of timesample,physical points, subslopes - REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol] - REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa] - REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa] - REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration -! OUTPUT - REAL,INTENT(OUT) :: Tcond(ngrid) ! Physical points : condensation temperature of CO2, yearly averaged - -! LOCAL - - INTEGER :: ig,it ! for loop - REAL :: ave ! intermediate to compute average - -!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - DO ig = 1,ngrid - ave = 0 - DO it = 1,timelen - ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(ig,it)*ps_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100)) - ENDDO - Tcond(ig) = ave/timelen - ENDDO -RETURN - - -end subroutine -end module Index: trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90 (revision 2993) +++ trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90 (revision 2995) @@ -25,5 +25,5 @@ USE comsoil_h_pem, only: soil_pem,fluxgeo,water_reservoir_nom,depth_breccia,depth_bedrock,reg_thprop_dependp USE adsorption_mod,only: adsorption_pem - USE co2glaciers_mod, only: co2glaciersflow + USE glaciers_mod, only: co2glaciersflow,h2oglaciersflow use ice_table_mod, only: icetable_equilibrium, icetable_dynamic @@ -79,4 +79,7 @@ co2glaciersflow = .true. CALL getin('co2glaciersflow', co2glaciersflow) + + h2oglaciersflow = .true. + CALL getin('h2oglaciersflow', h2oglaciersflow) reg_thprop_dependp = .false. Index: trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90 (revision 2995) +++ trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90 (revision 2995) @@ -0,0 +1,317 @@ +module glaciers_mod + + implicit none + LOGICAL co2glaciersflow ! True by default, to compute co2 ice flow. Read in pem.def + LOGICAL h2oglaciersflow ! True by default, to compute co2 ice flow. Read in pem.def +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Compute CO2 glacier flows +!!! +!!! Author: LL +!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +contains + + +subroutine co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do +!!! the ice transfer +!!! +!!! +!!! Author: LL +!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + +IMPLICIT NONE + +! arguments +! --------- + +! Inputs: + INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! number of time sample, physical points, subslopes, index of the flat subslope + REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles + REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol] + REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical x Time field: surface pressure given by the GCM [Pa] + REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surface pressure from the GCM [Pa] + REAL,INTENT(IN) :: global_ave_ps_PEM ! global averaged surface pressure during the PEM iteration [Pa] + +! Ouputs: + REAL,INTENT(INOUT) :: co2ice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] + REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes + REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid) ! same but within the mesh + + +! Local + REAL :: Tcond(ngrid,nslope) ! Physical field: CO2 condensation temperature [K] + REAL :: hmax(ngrid,nslope) ! Physical x Slope field: maximum thickness for co2 glacier before flow + +!----------------------------- + call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond) + + call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax) + + call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh) + RETURN +end subroutine + + + + + +subroutine h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do +!!! the ice transfer +!!! +!!! +!!! Author: LL +!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + +IMPLICIT NONE + +! arguments +! --------- + +! Inputs: + INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! number of time sample, physical points, subslopes, index of the flat subslope + REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles + REAL,INTENT(IN) :: Tice(ngrid,nslope) ! Ice Temperature [K] +! Ouputs: + REAL,INTENT(INOUT) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] + REAL,INTENT(INOUT) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes + REAL,INTENT(INOUT) :: flag_h2oflow_mesh(ngrid) ! same but within the mesh +! Local + REAL :: hmax(ngrid,nslope) ! Physical x Slope field: maximum thickness for co2 glacier before flow + +!----------------------------- + + call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax) + + call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh) + + RETURN +end subroutine + + + +subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax) + + USE comconst_mod, ONLY: pi,g + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle +!!! before initating flow +!!! +!!! Author: LL,based on work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022) +!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + IMPLICIT NONE + +! arguments +! -------- + +! Inputs + INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes + INTEGER,INTENT(IN) :: iflat ! index of the flat subslope + REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg] + REAL,INTENT(IN) :: Tice(ngrid,nslope) ! Physical field: ice temperature [K] + character(len=30), INTENT(IN) :: name_ice ! Nature of the ice +! Outputs + REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum thickness before flaw [m] +! Local + DOUBLE PRECISION :: tau_d ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith + REAL :: rho(ngrid,nslope) ! co2 ice density [kg/m^3] + INTEGER :: ig,islope ! loop variables + REAL :: slo_angle + + +! 1. Compute rho + if(name_ice.eq."co2") then + DO ig = 1,ngrid + DO islope = 1,nslope + rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3 ! Mangan et al. 2017 + tau_d = 5.e3 + ENDDO + ENDDO + elseif (name_ice.eq."h2o") then + DO ig = 1,ngrid + DO islope = 1,nslope + rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+ 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 + tau_d = 1.e5 + ENDDO + ENDDO + else + call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1) + endif + +! 3. Compute max thickness + DO ig = 1,ngrid + DO islope = 1,nslope + if(islope.eq.iflat) then + hmax(ig,islope) = 1.e8 + else + slo_angle = abs(def_slope_mean(islope)*pi/180.) + hmax(ig,islope) = tau_d/(rho(ig,islope)*g*slo_angle) + endif + ENDDO + ENDDO +RETURN + +end subroutine + + + + +subroutine transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Transfer the excess of ice from one subslope to another +!!! No transfer between mesh at the time +!!! Author: LL +!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + USE comconst_mod, ONLY: pi + + +implicit none + +! arguments +! -------- + +! Inputs + INTEGER, INTENT(IN) :: ngrid,nslope !# of physical points and subslope + INTEGER, INTENT(IN) :: iflat ! index of the flat subslope + REAL, INTENT(IN) :: subslope_dist(ngrid,nslope) ! Distribution of the subgrid slopes within the mesh + REAL, INTENT(IN) :: def_slope_mean(nslope) ! values of the subgrid slopes + REAL, INTENT(IN) :: hmax(ngrid,nslope) ! maximum height of the glaciers before initiating flow [m] + REAL, INTENT(IN) :: Tice(ngrid,nslope) ! Ice temperature[K] + character(len=30), INTENT(IN) :: name_ice ! Nature of the ice + +! Outputs + REAL, INTENT(INOUT) :: qice(ngrid,nslope) ! CO2 in the subslope [kg/m^2] + REAL, INTENT(INOUT) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope + REAL, INTENT(INOUT) :: flag_flowmesh(ngrid) ! boolean to check if there is flow in the mesh +! Local + INTEGER ig,islope ! loop + REAL rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3] + INTEGER iaval ! ice will be transfered here + + + +! 0. Compute rho + if(name_ice.eq."co2") then + DO ig = 1,ngrid + DO islope = 1,nslope + rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3 ! Mangan et al. 2017 + ENDDO + ENDDO + elseif (name_ice.eq."h2o") then + DO ig = 1,ngrid + DO islope = 1,nslope + rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+ 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 + ENDDO + ENDDO + else + call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1) + endif + +! 1. Compute the transfer of ice + + DO ig = 1,ngrid + DO islope = 1,nslope + IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground +! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely on flat ground + IF(qice(ig,islope).ge.rho(ig,islope)*hmax(ig,islope) * & + cos(pi*def_slope_mean(islope)/180.)) THEN +! Second: determine the flatest slopes possible: + IF(islope.gt.iflat) THEN + iaval=islope-1 + ELSE + iaval=islope+1 + ENDIF + do while ((iaval.ne.iflat).and. & + (subslope_dist(ig,iaval).eq.0)) + IF(iaval.gt.iflat) THEN + iaval=iaval-1 + ELSE + iaval=iaval+1 + ENDIF + enddo + qice(ig,iaval) = qice(ig,iaval) + & + (qice(ig,islope) - rho(ig,islope)*hmax(ig,islope) * & + cos(pi*def_slope_mean(islope)/180.)) * & + subslope_dist(ig,islope)/subslope_dist(ig,iaval) * & + cos(pi*def_slope_mean(iaval)/180.) / & + cos(pi*def_slope_mean(islope)/180.) + + qice(ig,islope)=rho(ig,islope)*hmax(ig,islope) * & + cos(pi*def_slope_mean(islope)/180.) + + flag_flow(ig,islope) = 1. + flag_flowmesh(ig) = 1. + ENDIF ! co2ice > hmax + ENDIF ! iflat + ENDDO !islope + ENDDO !ig +RETURN +end subroutine + + subroutine computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Compute CO2 condensation temperature +!!! +!!! Author: LL +!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2 + +implicit none + +! arguments: +! ---------- + +! INPUT + INTEGER,INTENT(IN) :: timelen, ngrid,nslope ! # of timesample,physical points, subslopes + REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol] + REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa] + REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa] + REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration +! OUTPUT + REAL,INTENT(OUT) :: Tcond(ngrid,nslope) ! Physical points : condensation temperature of CO2, yearly averaged + +! LOCAL + + INTEGER :: ig,it,islope ! for loop + REAL :: ave ! intermediate to compute average + +!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + DO ig = 1,ngrid + ave = 0 + DO it = 1,timelen + ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(ig,it)*ps_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100)) + ENDDO + DO islope = 1,nslope + Tcond(ig,islope) = ave/timelen + ENDDO + ENDDO +RETURN + + +end subroutine +end module Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 2993) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 2995) @@ -93,5 +93,5 @@ use conf_pem_mod, only: conf_pem use pemredem, only: pemdem0,pemdem1 - use co2glaciers_mod, only: co2glaciers_evol,co2glaciersflow + use glaciers_mod, only: co2glaciers_evol,h2oglaciersflow,co2glaciersflow,h2oglaciersflow 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, & @@ -942,4 +942,5 @@ call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) + call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) ENDDO