Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 3095) +++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 3098) @@ -37,5 +37,6 @@ use aeropacity_mod, only: iddist, topdustref USE mod_phys_lmdz_transfert_para, ONLY: bcast - USE paleoclimate_mod,ONLY: paleoclimate,albedo_perenialco2 + USE paleoclimate_mod,ONLY: paleoclimate,albedo_perenialco2, + & lag_layer use microphys_h, only: mteta @@ -338,8 +339,14 @@ ! PALEOCLIMATE - write(*,*)"Is it paleoclimate run?" + write(*,*)"Using lag layer??" + lag_layer=.false. + call getin_p("lag_layer",lag_layer) + write(*,*) " lag_layer = ", lag_layer + + write(*,*)"Is it paleoclimate run?" paleoclimate=.false. ! default value call getin_p("paleoclimate",paleoclimate) write(*,*)" paleoclimate = ",paleoclimate + write(*,*)"Albedo for perenial CO2 ice?" Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 (revision 3095) +++ trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 (revision 3098) @@ -477,4 +477,5 @@ call getin("subsurface_ice_depth",ice_depth) + z0(1) = z0_default ! default value for roughness write(*,*) 'Surface roughness length z0 (m)?' @@ -716,4 +717,8 @@ write(*,*) 'Prescribed atmospheric water vapor profile' write(*,*) 'Unless it reaches saturation (maximal value)' + else if (atm_wat_profile .eq. 2) then + write(*,*) 'Prescribed atmospheric water vapor profile' + write(*,*) 'like present day northern midlatitude' + write(*,*) 'Unless it reaches saturation (maximal value)' else error stop 'Water vapor profile value not correct!' Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 (revision 3095) +++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 (revision 3098) @@ -179,5 +179,17 @@ endif endif - endif + !EV profile +! IF(atm_wat_profile.eq.2) THEN +! DO ilayer=1,16 +! q(1,ilayer,igcm_h2o_vap)=min(zqsat(ilayer),6e-5) +! ENDDO! ilayer=1,16 +! DO ilayer=17,22 +! q(1,ilayer,igcm_h2o_vap)=min(zqsat(ilayer),3e-5) +! ENDDO! ilayer=17,22 +! DO ilayer=23,nlayer +! q(1,ilayer,igcm_h2o_vap)=0 +! ENDDO! ilayer=23,nlayer +! endif + endif ! Call physics Index: trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90 (revision 3095) +++ trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90 (revision 3098) @@ -15,8 +15,10 @@ !$OMP THREADPRIVATE(paleoclimate) - real, save, allocatable :: lag_h2o_ice(:,:) ! Thickness of the lag before H2O ice [m] + real, save, allocatable :: h2o_ice_depth(:,:) ! Thickness of the lag before H2O ice [m] real, save, allocatable :: lag_co2_ice(:,:) ! Thickness of the lag before CO2 ice [m] real, save :: albedo_perenialco2 ! Albedo for perenial co2 ice [1] -!$OMP THREADPRIVATE(lag_h2o_ice,lag_co2_ice,albedo_perenialco2) + real, save, allocatable :: d_coef(:,:) ! Diffusion coeficent + LOGICAL,SAVE :: lag_layer ! does lag layer is present? +!$OMP THREADPRIVATE(h2o_ice_depth,lag_co2_ice,albedo_perenialco2) CONTAINS @@ -29,6 +31,7 @@ integer,intent(in) :: nslope ! number of slope within a mesh - allocate(lag_h2o_ice(ngrid,nslope)) + allocate(h2o_ice_depth(ngrid,nslope)) allocate(lag_co2_ice(ngrid,nslope)) + allocate(d_coef(ngrid,nslope)) end subroutine ini_paleoclimate_h @@ -36,5 +39,6 @@ implicit none - if (allocated(lag_h2o_ice)) deallocate(lag_h2o_ice) + if (allocated(d_coef)) deallocate(d_coef) + if (allocated(h2o_ice_depth)) deallocate(h2o_ice_depth) if (allocated(lag_co2_ice)) deallocate(lag_co2_ice) end subroutine end_paleoclimate_h Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90 (revision 3095) +++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90 (revision 3098) @@ -27,5 +27,5 @@ use comsoil_h, only: flux_geo USE comslope_mod, ONLY: nslope, major_slope - USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice + USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef USE comcstfi_h, only: pi use geometry_mod, only: latitude @@ -744,10 +744,19 @@ if (startphy_file) then ! Depth of H2O lag - call get_field("lag_h2o_ice",lag_h2o_ice,found,indextime) - if (.not.found) then - write(*,*) "phyetat0: Failed loading : ", & - " is set as -1 (no subsurface ice)" - lag_h2o_ice(:,:) = -1. - endif + call get_field("h2o_ice_depth",h2o_ice_depth,found,indextime) + if (.not.found) then + write(*,*) "phyetat0: Failed loading : ", & + " is set as -1 (no subsurface ice)" + h2o_ice_depth(:,:) = -1. + endif + +! Diffuesion coeficent + call get_field("d_coef",d_coef,found,indextime) + if (.not.found) then + write(*,*) "phyetat0: Failed loading : ", & + " is set as 4e-4 (defualt value)" + d_coef(:,:) = 4e-4 + endif + ! Depth of CO2 lag @@ -774,7 +783,8 @@ endif ! notfound else ! no startfiphyle - lag_h2o_ice(:,:) = -1. + h2o_ice_depth(:,:) = -1. lag_co2_ice(:,:) = -1. perenial_co2ice(:,:) = 0. + d_coef(:,:)= 4.e-4 if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then DO islope = 1,nslope @@ -786,7 +796,8 @@ endif !startphy_file else - lag_h2o_ice(:,:) = -1. + h2o_ice_depth(:,:) = -1. lag_co2_ice(:,:) = -1. perenial_co2ice(:,:) = 0. + d_coef(:,:)= 4.e-4 endif !paleoclimate Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/phyredem.F90 (revision 3095) +++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90 (revision 3098) @@ -25,5 +25,5 @@ use time_phylmdz_mod, only: daysec use comslope_mod, ONLY: nslope - USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice + USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice implicit none @@ -156,5 +156,5 @@ ! Paleoclimate outputs if(paleoclimate) then - call put_field("lag_h2o_ice","Depth of the H2O lags",lag_h2o_ice) + call put_field("h2o_ice_depth","Depth to the fisrt H2O ice",h2o_ice_depth) call put_field("lag_co2_ice","Depth of the CO2 lags",lag_co2_ice) endif Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 3095) +++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 3098) @@ -1514,5 +1514,6 @@ & zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th, & zcondicea_co2microp,sensibFlux, - & dustliftday,local_time,watercap,dwatercap_dif) + & dustliftday,local_time,watercap,dwatercap_dif, + & tsoil,nsoilmx) DO ig=1,ngrid Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F (revision 3095) +++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F (revision 3098) @@ -13,6 +13,6 @@ $ pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true, $ hfmax,pcondicea_co2microp,sensibFlux, - $ dustliftday,local_time,watercap, dwatercap_dif) - + $ dustliftday,local_time,watercap, dwatercap_dif, + $ tsoil,nsoil) use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number, & igcm_dust_submicron, igcm_h2o_vap, @@ -32,4 +32,6 @@ & subslope_dist,major_slope,iflat use microphys_h, only: To + use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer + use comsoil_h, only: layer, mlayer IMPLICIT NONE @@ -61,5 +63,5 @@ c ---------- - INTEGER,INTENT(IN) :: ngrid,nlay + INTEGER,INTENT(IN) :: ngrid,nlay,nsoil REAL,INTENT(IN) :: ptimestep REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1) @@ -75,4 +77,5 @@ REAL,INTENT(IN) :: pcapcal(ngrid,nslope) REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1) + REAL,INTENT(IN) :: tsoil(ngrid,nsoil,nslope) c Argument added for condensation: @@ -102,5 +105,5 @@ REAL :: pt(ngrid,nlay) - INTEGER ilev,ig,ilay,nlev,islope + INTEGER ilev,ig,ilay,nlev,islope,ik,lice REAL z4st,zdplanck(ngrid) @@ -118,4 +121,11 @@ REAL zcst1 REAL zu2(ngrid) + REAL Tice(ngrid,nslope) ! subsurface temperature where ice is located. + REAL qeq(ngrid,nslope) ! saturation water vapor in the subsurface + REAL dist_up(ngrid,nslope) !distance from ice to layer above + REAL dist_down(ngrid,nslope) !distance from ice to layer down + REAL dist_sum(ngrid,nslope) ! sum of distance + REAL zdqsdif_ssi(ngrid,nslope) !SSI flux + REAL zdqsdif_ssi_frost(ngrid,nslope) !SSI-frost interaction EXTERNAL SSUM,SCOPY @@ -155,4 +165,6 @@ REAL rho(ngrid) ! near surface air density REAL qsat(ngrid) + REAL qsat2(ngrid,nslope) + REAL resist(ngrid,nslope) !subsurface ice flux reduction coef REAL hdoflux(ngrid,nslope) ! value of vapour flux of HDO @@ -272,5 +284,5 @@ zdqsdif(1:ngrid)=0 dwatercap_dif(1:ngrid,1:nslope)=0 - +! h2o_ice_depth(1:ngrid,1:nslope)=5 c ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa @@ -914,5 +926,5 @@ c Utilisation d'un sous pas de temps afin c de decrire le flux de chaleur latente - + do iq=1,nq if ((water).and.(iq.eq.igcm_h2o_vap)) then @@ -920,4 +932,7 @@ DO islope = 1,nslope DO ig=1,ngrid + + + zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/ & cos(pi*def_slope_mean(islope)/180.) @@ -943,5 +958,4 @@ zq_tmp_vap(ig,:,:) =zq(ig,:,:) c Debut du sous pas de temps - DO tsub=1,nsubtimestep(ig) @@ -980,7 +994,89 @@ if ((-zdqsdif(ig)*subtimestep) & .gt.(zqsurf(ig))) then + + + !EV subsurface ice + IF(h2o_ice_depth(ig,islope) .gt. 0 .and. lag_layer) + & then + zdqsdif(ig)= + & -zqsurf(ig)/subtimestep + zqsurf(ig)=0 + + DO ik=0,nsoil-2 ! go through all the layers to find the ice locations + IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and. + & (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN + lice = ik+1 + EXIT + ENDIF + ENDDO !of subsurface loop + IF (lice .gt. 1) then !calculate the distance from the layers + dist_up(ig,islope)=(h2o_ice_depth(ig,islope) + & -mlayer(lice-1)) + dist_down(ig,islope)=(mlayer(lice) + & -h2o_ice_depth(ig,islope)) + dist_sum(ig,islope)=dist_up(ig,islope) + & +dist_down(ig,islope) + Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp + & *tsoil(ig,lice-1,islope) + & /dist_sum(ig,islope))+ + & (dist_down(ig,islope)*tsoil(ig,lice,islope) + & /dist_sum(ig,islope)) + ELSE + Tice(ig,islope)=tsoil(ig,1,islope) + ENDIF + call watersat(1,Tice(ig,1),pplev(ig,1) + & ,qsat2(ig,1)) + qeq(ig,1)=(ztsrf(ig)/Tice(ig,1)) + & *qsat2(ig,1) + resist(ig,1)=(1+(h2o_ice_depth(ig,islope)*zcdh(ig) + & /d_coef(ig,islope))) + !write(*,*)'R=',resist(ig,islope) + !write(*,*)'zice=',h2o_ice_depth(ig,islope) + !write(*,*)'D=',d_coef(ig,islope) + !write(*,*)'zcdh=',zcdh(ig) + zb(ig,1)=zb(ig,1)/resist(ig,1) ! change zb to account subsurface ice + !vdifc algorithem !!!!!needs to change reseist io to the mean + !beacuse of the slopes!!!!! + z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) + zc(ig,nlay)=za(ig,nlay)*zq_tmp_vap(ig,nlay,iq)*z1(ig) + zd(ig,nlay)=zb(ig,nlay)*z1(ig) + DO ilay=nlay-1,2,-1 + z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ + $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) + zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+ + $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) + zd(ig,ilay)=zb(ig,ilay)*z1(ig) + ENDDO + z1(ig)=1./(za(ig,1)+zb(ig,1)+ + $ zb(ig,2)*(1.-zd(ig,2))) + zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+ + $ zb(ig,2)*zc(ig,2)) * z1(ig) + zd(ig,1)=zb(ig,1)*z1(ig) + zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig) + zdqsdif_ssi(ig,1)=rho(ig)*dryness(ig)*zcdv(ig) + & *(zq1temp(ig)-qeq(ig,islope)) + call write_output('zdq_ssi', + & '','',zdqsdif_ssi(ig,1)) + call write_output('qeq', + & '','',qeq(ig,1)) + call write_output('q1', + & '','',zq1temp(ig)) + !write(*,*)'qeq=',qeq(ig,1) + !write(*,*)'q1=',zq1temp(ig) + !write(*,*)'zdqsdif_ssi=',zdqsdif_ssi(ig,1) + !I should some all the interactions with the SSI accoridng + !to the statistics and the evaluate + + !write(*,*) "zdq_ssi*t", zdqsdif_ssi(ig)*subtimestep + if (zdqsdif_ssi(ig,1)<0) then + !zdqsdif(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep) + zdqsdif(ig)=zdqsdif(ig)+zdqsdif_ssi(ig,1) + call write_output('zdq_zdqssi', + & '','',zdqsdif(ig)+zdqsdif_ssi(ig,1)) + endif + ELSE c pdqsdif > 0 : ice condensing c pdqsdif < 0 : ice subliming -c write(*,*) "subliming more than available frost: qsurf!" + write(*,*) "subliming more than available frost: qsurf!" zdqsdif(ig)= & -zqsurf(ig)/subtimestep @@ -991,11 +1087,37 @@ $ (-zdqsdif(ig)) *subtimestep) *z1(ig) zq1temp(ig)=zc(ig,1) + ENDIF !if h2o_ice_depth>0 and lag_layer endif !if .not.watercaptag(ig) endif ! if sublim more than surface + + ! subsurface ice <--> frost interaction !EV + IF(h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer + & .and. zqsurf(ig) .gt. 0) then + zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1) + & /h2o_ice_depth(ig,1)) + & *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1)) + !needs to change to the mean of eq + zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1) + !!!!! zdsqdif_ssi_frosst need to be changed to an + !average + ELSEIF (h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer + & .and. watercaptag(ig)) then + zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1) + & /h2o_ice_depth(ig,1)) + & *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1)) + zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1) + !needs to change to the mean of eq + ENDIF + call write_output('zdqsdif_ssi_frost', + & '','',zdqsdif_ssi_frost(ig,1)) + call write_output('subtimestep', + & '','',subtimestep) + ! ENDDO !subsurface ice subslope + + c Starting upward calculations for water : c Actualisation de h2o_vap dans le premier niveau zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig) - c Take into account the H2O latent heat impact on the surface temperature if (latentheat_surfwater) then @@ -1009,5 +1131,4 @@ & *zq_tmp_vap(ig,ilay-1,iq) ENDDO - c Subtimestep water budget : @@ -1016,5 +1137,7 @@ zqsurf(ig)= zqsurf(ig)+( & zdqsdif(ig))*subtimestep - + if (zqsurf(ig)<0) then + zqsurf(ig)=0 + endif c Monitoring instantaneous latent heat flux in W.m-2 : zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+ @@ -1033,5 +1156,4 @@ c Fin du sous pas de temps ENDDO ! tsub=1,nsubtimestep - c Integration of subtimestep temp and water budget : c (btw could also compute the post timestep temp and ice @@ -1044,5 +1166,4 @@ & cos(pi*def_slope_mean(islope)/180.)) & /ptimestep - c if subliming more than qsurf(ice) and on watercaptag, water c sublimates from watercap reservoir (dwatercap_dif is <0) @@ -1063,5 +1184,4 @@ ENDDO ! of DO ig=1,ngrid ENDDO ! islope - c Some grid box averages: interface with the atmosphere DO ig = 1,ngrid @@ -1075,5 +1195,4 @@ ENDDO ENDDO - ! Recompute values in kg/m^2 slopped DO ig = 1,ngrid @@ -1167,5 +1286,10 @@ & "Ground ice latent heat flux", & "W.m-2",surf_h2o_lh(:,iflat)) - + call write_output('zdq_subtimestep', + & '','',zdqsdif(:)*subtimestep) + call write_output('zdq_end', + & '','',zdqsdif(:)) + call write_output('vdifc_subtimestep', + & '','',subtimestep) C Diagnostic output for HDO ! if (hdo) then