Index: trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F (revision 2952) +++ trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F (revision 2953) @@ -8,5 +8,5 @@ CONTAINS - SUBROUTINE co2condens(ngrid,nlayer,nq,ptimestep, + SUBROUTINE co2condens(ngrid,nlayer,nq,nslope,ptimestep, $ pcapcal,pplay,pplev,ptsrf,pt, $ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, @@ -34,4 +34,5 @@ USE vertical_layers_mod, ONLY: ap, bp #endif + use comslope_mod, ONLY: subslope_dist,def_slope_mean IMPLICIT NONE c======================================================================= @@ -59,10 +60,11 @@ INTEGER,INTENT(IN) :: nlayer ! number of vertical layers INTEGER,INTENT(IN) :: nq ! number of tracers + INTEGER,INTENT(IN) :: nslope ! number of subslope REAL,INTENT(IN) :: ptimestep ! physics timestep (s) - REAL,INTENT(IN) :: pcapcal(ngrid) + REAL,INTENT(IN) :: pcapcal(ngrid,nslope) REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa) REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) - REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K) + REAL,INTENT(IN) :: ptsrf(ngrid,nslope) ! surface temperature (K) REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K) REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2) @@ -73,5 +75,5 @@ REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2) ! from previous physical processes - REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from + REAL,INTENT(IN) :: pdtsrf(ngrid,nslope) ! tendency on surface temperature from ! previous physical processes (K/s) REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s) @@ -84,12 +86,12 @@ REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1) - REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2) - REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface - REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface + REAL,INTENT(INOUT) :: piceco2(ngrid,nslope) ! CO2 ice on the surface (kg.m-2) + REAL,INTENT(INOUT) :: psolaralb(ngrid,2,nslope) ! albedo of the surface + REAL,INTENT(INOUT) :: pemisurf(ngrid,nslope) ! emissivity of the surface REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius ! tendencies due to CO2 condensation/sublimation: REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s) - REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s) + REAL,INTENT(OUT) :: pdtsrfc(ngrid,nslope) ! tendency on surface temperature (K/s) REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s) REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2) @@ -112,7 +114,9 @@ REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm) REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface) - REAL zdiceco2(ngrid) + REAL zdiceco2(ngrid,nslope) + REAL zdiceco2_mesh_avg(ngrid) REAL zcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s) - REAL zcondices(ngrid) ! condensation rate on the ground (kg/m2/s) + REAL zcondices(ngrid,nslope) ! condensation rate on the ground (kg/m2/s) + REAL zcondices_mesh_avg(ngrid) ! condensation rate on the ground (kg/m2/s) REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s) REAL condens_layer(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s) @@ -122,5 +126,5 @@ REAL zu(nlayer),zv(nlayer) REAL zqc(nlayer,nq),zq1(nlayer) - REAL ztsrf(ngrid) + REAL ztsrf(ngrid,nslope) REAL ztc(nlayer), ztm(nlayer+1) REAL zum(nlayer+1) , zvm(nlayer+1) @@ -129,7 +133,7 @@ REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix REAL availco2 - LOGICAL condsub(ngrid) - - real :: emisref(ngrid) + LOGICAL condsub(ngrid,nslope) + + real :: emisref(ngrid,nslope) REAL zdq_scav(ngrid,nlayer,nq) ! tendency due to scavenging by co2 @@ -170,4 +174,16 @@ REAL masseq(nlayer),wq(nlayer+1) INTEGER ifils,iq2 + +c Subslope: + + REAL :: emisref_tmp(ngrid) + REAL :: alb_tmp(ngrid,2) ! local + REAL :: zcondices_tmp(ngrid) ! local + REAL :: piceco2_tmp(ngrid) ! local + REAL :: pemisurf_tmp(ngrid)! local + LOGICAL :: condsub_tmp(ngrid) !local + REAL :: zfallice_tmp(ngrid,nlayer+1) + REAL :: condens_layer_tmp(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s) + INTEGER :: islope c---------------------------------------------------------------------- @@ -231,9 +247,11 @@ pdqc(1:ngrid,1:nlayer,1:nq) = 0 - zcondices(1:ngrid) = 0. - pdtsrfc(1:ngrid) = 0. + zcondices(1:ngrid,1:nslope) = 0. + zcondices_mesh_avg(1:ngrid)=0. + pdtsrfc(1:ngrid,1:nslope) = 0. pdpsrf(1:ngrid) = 0. - condsub(1:ngrid) = .false. - zdiceco2(1:ngrid) = 0. + condsub(1:ngrid,1:nslope) = .false. + zdiceco2(1:ngrid,1:nslope) = 0. + zdiceco2_mesh_avg(1:ngrid)=0. zfallheat=0 @@ -301,5 +319,5 @@ pdtc(ig,l)=0. IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN - condsub(ig)=.true. + condsub(ig,:)=.true. IF (zfallice(ig,l+1).gt.0) then zfallheat=zfallice(ig,l+1)* @@ -347,5 +365,9 @@ condens_column(ig) = sum(condens_layer(ig,:)) zfallice(ig,1) = zdqssed_co2(ig) - piceco2(ig) = piceco2(ig) + zdqssed_co2(ig)*ptimestep + DO islope = 1,nslope + piceco2(ig,islope) = piceco2(ig,islope) + + & zdqssed_co2(ig)*ptimestep * + & cos(pi*def_slope_mean(islope)/180.) + ENDDO ENDDO call write_output("co2condens_zfallice", @@ -371,5 +393,8 @@ ztcondsol(ig)= & 1./(bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1))) - ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep + DO islope = 1,nslope + ztsrf(ig,islope) = ptsrf(ig,islope) + + & pdtsrf(ig,islope)*ptimestep + ENDDO ENDDO @@ -388,11 +413,19 @@ icap=1 ENDIF + + DO islope = 1,nslope +c Need first to put piceco2_slope(ig,islope) in kg/m^2 flat + + piceco2(ig,islope) = piceco2(ig,islope) + & /cos(pi*def_slope_mean(islope)/180.) + c c Loop on where we have condensation/ sublimation - IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond + IF ((ztsrf(ig,islope) .LT. ztcondsol(ig)) .OR. ! ground cond $ (zfallice(ig,1).NE.0.) .OR. ! falling snow - $ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim. - $ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN - condsub(ig) = .true. + $ ((ztsrf(ig,islope) .GT. ztcondsol(ig)) .AND. ! ground sublim. + $ ((piceco2(ig,islope)+zfallice(ig,1)*ptimestep) + $ .NE. 0.))) THEN + condsub(ig,islope) = .true. IF (zfallice(ig,1).gt.0 .and. .not. co2clouds) then @@ -406,7 +439,20 @@ c condensation or partial sublimation of CO2 ice c """"""""""""""""""""""""""""""""""""""""""""""" - zcondices(ig) = pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) - & /(latcond*ptimestep) - zfallheat - pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep + if(ztsrf(ig,islope).LT. ztcondsol(ig)) then +c condensation + zcondices(ig,islope)=pcapcal(ig,islope) + & *(ztcondsol(ig)-ztsrf(ig,islope)) +c & /((latcond+cpp*(zt(ig,1)-ztcondsol(ig)))*ptimestep) + & /(latcond*ptimestep) + & - zfallheat + else +c sublimation + zcondices(ig,islope)=pcapcal(ig,islope) + & *(ztcondsol(ig)-ztsrf(ig,islope)) + & /(latcond*ptimestep) + & - zfallheat + endif + pdtsrfc(ig,islope) = (ztcondsol(ig) - ztsrf(ig,islope)) + & /ptimestep #ifdef MESOSCALE print*, "not enough CO2 tracer in 1st layer to condense" @@ -421,11 +467,13 @@ availco2 = pq(ig,1,ico2)*((ap(1)-ap(2))+ & (bp(1)-bp(2))*(pplev(ig,1)/g - - & (zcondices(ig) + zfallice(ig,1))*ptimestep)) - if ((zcondices(ig) + condens_layer(ig,1))*ptimestep + & (zcondices(ig,islope) + zfallice(ig,1)) + & *ptimestep)) + if ((zcondices(ig,islope) + condens_layer(ig,1)) + & *ptimestep & .gt.availco2) then - zcondices(ig) = availco2/ptimestep - + zcondices(ig,islope) = availco2/ptimestep - & condens_layer(ig,1) - pdtsrfc(ig) = (latcond/pcapcal(ig))* - & (zcondices(ig)+zfallheat) + pdtsrfc(ig,islope) = (latcond/pcapcal(ig,islope))* + & (zcondices(ig,islope)+zfallheat) end if #else @@ -437,35 +485,52 @@ #endif -c If the entire CO2 ice layer sublimes +c If the entire CO2 ice layer sublimes on the slope c """""""""""""""""""""""""""""""""""""""""""""""""""" c (including what has just condensed in the atmosphere) if (co2clouds) then - IF((piceco2(ig)/ptimestep).LE. - & -zcondices(ig))THEN - zcondices(ig) = -piceco2(ig)/ptimestep - pdtsrfc(ig)=(latcond/pcapcal(ig)) * - & (zcondices(ig)+zfallheat) + IF((piceco2(ig,islope)/ptimestep).LE. + & -zcondices(ig,islope))THEN + zcondices(ig,islope) = -piceco2(ig,islope)/ptimestep + pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) * + & (zcondices(ig,islope)+zfallheat) END IF else - IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE. - & -zcondices(ig))THEN - zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) - pdtsrfc(ig)=(latcond/pcapcal(ig)) * - & (zcondices(ig)+zfallheat) + IF((piceco2(ig,islope)/ptimestep+zfallice(ig,1)).LE. + & -zcondices(ig,islope))THEN + zcondices(ig,islope) = -piceco2(ig,islope) + & /ptimestep - zfallice(ig,1) + pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) * + & (zcondices(ig,islope)+zfallheat) END IF end if -c Changing CO2 ice amount and pressure : +c Changing CO2 ice amount and pressure per slope: c """""""""""""""""""""""""""""""""""" - zdiceco2(ig) = zcondices(ig) + zfallice(ig,1) + zdiceco2(ig,islope) = zcondices(ig,islope)+zfallice(ig,1) & + condens_column(ig) if (co2clouds) then ! add here only direct condensation/sublimation - piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep + piceco2(ig,islope) = piceco2(ig,islope) + + & zcondices(ig,islope)*ptimestep else ! add here also CO2 ice in the atmosphere - piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep + piceco2(ig,islope) = piceco2(ig,islope) + + & zdiceco2(ig,islope)*ptimestep end if - pdpsrf(ig) = -zdiceco2(ig)*g + + zcondices_mesh_avg(ig) = zcondices_mesh_avg(ig) + + & zcondices(ig,islope)* subslope_dist(ig,islope) + + zdiceco2_mesh_avg(ig) = zdiceco2_mesh_avg(ig) + + & zdiceco2(ig,islope)* subslope_dist(ig,islope) + + END IF ! if there is condensation/sublimation + + piceco2(ig,islope) = piceco2(ig,islope) + & * cos(pi*def_slope_mean(islope)/180.) + + ENDDO !islope + + pdpsrf(ig) = -zdiceco2_mesh_avg(ig)*g IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN @@ -480,6 +545,7 @@ & 'condensing more than total mass', 1) ENDIF - END IF ! if there is condensation/sublimation + ENDDO ! of DO ig=1,ngrid + c ******************************************************************** @@ -489,20 +555,32 @@ ! Check that amont of CO2 ice is not problematic DO ig=1,ngrid - if(.not.piceco2(ig).ge.0.) THEN - if(piceco2(ig).le.-5.e-8) print*, - $ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig) - piceco2(ig)=0. + DO islope = 1,nslope + if(.not.piceco2(ig,islope).ge.0.) THEN + if(piceco2(ig,islope).le.-5.e-8) print*, + $ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig,islope) + piceco2(ig,islope)=0. endif + ENDDO ENDDO ! Set albedo and emissivity of the surface ! ---------------------------------------- - CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref) + DO islope = 1,nslope + piceco2_tmp(:) = piceco2(:,islope) + alb_tmp(:,:) = psolaralb(:,:,islope) + emisref_tmp(:) = 0. + CALL albedocaps(zls,ngrid,piceco2_tmp,alb_tmp,emisref_tmp) + psolaralb(:,1,islope) = alb_tmp(:,1) + psolaralb(:,2,islope) = alb_tmp(:,2) + emisref(:,islope) = emisref_tmp(:) + ENDDO ! set pemisurf() to emissiv when there is bare surface (needed for co2snow) DO ig=1,ngrid - if (piceco2(ig).eq.0) then - pemisurf(ig)=emissiv - endif + DO islope = 1,nslope + if (piceco2(ig,islope).eq.0) then + pemisurf(ig,islope)=emissiv + endif + ENDDO ENDDO @@ -515,5 +593,5 @@ DO ig=1,ngrid - if (condsub(ig)) then + if (any(condsub(ig,:))) then do l=1,nlayer ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep @@ -527,5 +605,5 @@ c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up) c """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" - zmflux(1) = -zcondices(ig) - zdqssed_co2(ig) + zmflux(1) = -zcondices_mesh_avg(ig) - zdqssed_co2(ig) DO l=1,nlayer zmflux(l+1) = zmflux(l) - condens_layer(ig,l) @@ -702,6 +780,21 @@ c CO2 snow / clouds scheme c *************************************************************** - call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev, - & condens_layer,zcondices,zfallice,pemisurf) + DO islope = 1,nslope + emisref_tmp(:) = emisref(:,islope) + condsub_tmp(:) = condsub(:,islope) + condens_layer_tmp(:,:) = condens_layer(:,:)* + & cos(pi*def_slope_mean(islope)/180.) + zcondices_tmp(:) = zcondices(:,islope)* + & cos(pi*def_slope_mean(islope)/180.) + zfallice_tmp(:,:) = zfallice(:,:)* + & cos(pi*def_slope_mean(islope)/180.) + pemisurf_tmp(:) = pemisurf(:,islope) + + call co2snow(ngrid,nlayer,ptimestep,emisref_tmp,condsub_tmp, + & pplev,condens_layer_tmp,zcondices_tmp,zfallice_tmp, + & pemisurf_tmp) + pemisurf(:,islope) = pemisurf_tmp(:) + + ENDDO c *************************************************************** c Ecriture des diagnostiques @@ -739,5 +832,8 @@ & 1./(bcond-acond*log(.01*vmr_co2(ngrid,1)* & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) - pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep + DO islope = 1,nslope + pdtsrfc(ngrid,islope)=(ztcondsol(ngrid)- + & ztsrf(ngrid,islope))/ptimestep + ENDDO ! islope = 1,nslope endif endif Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 2952) +++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 2953) @@ -675,5 +675,4 @@ c initialize tracers c ~~~~~~~~~~~~~~~~~~ - CALL initracer(ngrid,nq,qsurf) @@ -681,7 +680,5 @@ c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CALL surfini(ngrid,qsurf) - CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) - c initialize soil c ~~~~~~~~~~~~~~~ @@ -1170,11 +1167,11 @@ ELSE ! not calling subslope, nslope might be > 1 DO islope = 1,nslope - IF(def_slope_mean(islope).ge.0.) THEN - sl_psi = 0. !Northward slope - ELSE - sl_psi = 180. !Southward slope - ENDIF sl_the=abs(def_slope_mean(islope)) IF (sl_the .gt. 1e-6) THEN + IF(def_slope_mean(islope).ge.0.) THEN + psi_sl(:) = 0. !Northward slope + ELSE + psi_sl(:) = 180. !Southward slope + ENDIF DO ig=1,ngrid ztim1=fluxsurf_dn_sw(ig,1,islope) @@ -1232,5 +1229,4 @@ ENDDO - c Net atmospheric radiative heating rate (K.s-1) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1254,4 +1250,5 @@ c Net radiative surface flux (W.m-2) c ~~~~~~~~~~~~~~~~~~~~~~~~~~ + c DO ig=1,ngrid @@ -1312,7 +1309,8 @@ c for radiative transfer & clearatm,icount,zday,zls, - & tsurf,qsurf_meshavg(:,igcm_co2),igout, - & totstormfract,tauscaling, + & tsurf_meshavg,qsurf_meshavg(:,igcm_co2), + & igout,totstormfract,tauscaling, & dust_rad_adjust,IRtoVIScoef, + & albedo_meshavg,emis_meshavg, c input sub-grid scale cloud & clearsky,totcloudfrac, @@ -1378,5 +1376,6 @@ & qsurf_meshavg(:,igcm_co2), & igout,aerosol,tauscaling,dust_rad_adjust, - & IRtoVIScoef,totstormfract,clearatm, + & IRtoVIScoef,albedo_meshavg,emis_meshavg, + & totstormfract,clearatm, & clearsky,totcloudfrac, & nohmons, @@ -1413,5 +1412,4 @@ endif ! end of if (dustinjection.gt.0) - c----------------------------------------------------------------------- c 4. Gravity wave and subgrid scale topography drag : @@ -1471,4 +1469,5 @@ ENDDO ENDIF + c ---------------------- @@ -2000,4 +1999,5 @@ c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep c Zdqssed isn't + call callsedim(ngrid,nlayer,ptimestep, & zplev,zzlev,zzlay,pt,pdt, @@ -2180,5 +2180,5 @@ zdqc(:,:,:) = 0. zdqsc(:,:,:) = 0. - CALL co2condens(ngrid,nlayer,nq,ptimestep, + CALL co2condens(ngrid,nlayer,nq,nslope,ptimestep, $ capcal,zplay,zplev,tsurf,pt, $ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, Index: trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90 (revision 2952) +++ trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90 (revision 2953) @@ -27,5 +27,5 @@ tsurf,co2ice,igout,totstormfract, & tauscaling,dust_rad_adjust, & - IRtoVIScoef, & + IRtoVIScoef,albedo,emis, & ! input sub-grid scale cloud clearsky,totcloudfrac, & @@ -40,7 +40,7 @@ ,rho_dust USE comcstfi_h, only: r,g,cpp,rcp - USE dimradmars_mod, only: albedo,naerkind + USE dimradmars_mod, only: naerkind USE comsaison_h, only: dist_sol,mu0,fract - USE surfdat_h, only: emis,zmea, zstd, zsig, hmons + USE surfdat_h, only: zmea, zstd, zsig, hmons USE callradite_mod, only: callradite use write_output_mod, only: write_output @@ -78,4 +78,6 @@ REAL, INTENT(IN) :: zls REAL, INTENT(IN) :: tsurf(ngrid) + REAL, INTENT(IN) :: albedo(ngrid,2) + REAL, INTENT(IN) :: emis(ngrid) REAL,INTENT(IN) :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2) INTEGER, INTENT(IN) :: igout @@ -255,4 +257,5 @@ ! 1. Call the second radiative transfer for stormdust, obtain the extra heating ! ********************************************************************* + CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, & emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, & Index: trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90 (revision 2952) +++ trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90 (revision 2953) @@ -22,5 +22,5 @@ icount,zday,zls,tsurf,co2ice,igout, & aerosol,tauscaling,dust_rad_adjust, & - IRtoVIScoef, & + IRtoVIScoef,albedo,emis, & ! input sub-grid scale rocket dust storm totstormfract,clearatm, & @@ -37,7 +37,7 @@ ,rho_dust USE comcstfi_h, only: r,g,cpp,rcp - USE dimradmars_mod, only: albedo,naerkind + USE dimradmars_mod, only: naerkind USE comsaison_h, only: dist_sol,mu0,fract - USE surfdat_h, only: emis,hmons,summit,alpha_hmons, & + USE surfdat_h, only: hmons,summit,alpha_hmons, & hsummit,contains_mons USE callradite_mod, only: callradite @@ -74,4 +74,6 @@ REAL, INTENT(IN) :: zls REAL, INTENT(IN) :: tsurf(ngrid) + REAL, INTENT(IN) :: albedo(ngrid,2) + REAL, INTENT(IN) :: emis(ngrid) REAL,INTENT(IN) :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2) INTEGER, INTENT(IN) :: igout Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F (revision 2952) +++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F (revision 2953) @@ -21,5 +21,5 @@ & igcm_stormdust_mass, igcm_stormdust_number use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness - USE comcstfi_h, ONLY: cpp, r, rcp, g + USE comcstfi_h, ONLY: cpp, r, rcp, g, pi use watersat_mod, only: watersat use turb_mod, only: turb_resolved, ustar, tstar @@ -29,4 +29,6 @@ use dust_param_mod, only: doubleq, submicron, lifting use write_output_mod, only: write_output + use comslope_mod, ONLY: nslope,def_slope,def_slope_mean, + & subslope_dist,major_slope,iflat IMPLICIT NONE @@ -65,11 +67,11 @@ REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) REAL,INTENT(IN) :: ph(ngrid,nlay) - REAL,INTENT(IN) :: ptsrf(ngrid),pemis(ngrid) + REAL,INTENT(IN) :: ptsrf(ngrid,nslope),pemis(ngrid,nslope) REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay) REAL,INTENT(IN) :: pdhfi(ngrid,nlay) - REAL,INTENT(IN) :: pfluxsrf(ngrid) + REAL,INTENT(IN) :: pfluxsrf(ngrid,nslope) REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay) - REAL,INTENT(OUT) :: pdtsrf(ngrid),pdhdif(ngrid,nlay) - REAL,INTENT(IN) :: pcapcal(ngrid) + REAL,INTENT(OUT) :: pdtsrf(ngrid,nslope),pdhdif(ngrid,nlay) + REAL,INTENT(IN) :: pcapcal(ngrid,nslope) REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1) @@ -87,9 +89,9 @@ c Traceurs : integer,intent(in) :: nq - REAL,INTENT(IN) :: pqsurf(ngrid,nq) + REAL,INTENT(IN) :: pqsurf(ngrid,nq,nslope) REAL :: zqsurf(ngrid) ! temporary water tracer real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) real,intent(out) :: pdqdif(ngrid,nlay,nq) - real,intent(out) :: pdqsdif(ngrid,nq) + real,intent(out) :: pdqsdif(ngrid,nq,nslope) REAL,INTENT(in) :: dustliftday(ngrid) REAL,INTENT(in) :: local_time(ngrid) @@ -100,5 +102,5 @@ REAL :: pt(ngrid,nlay) - INTEGER ilev,ig,ilay,nlev + INTEGER ilev,ig,ilay,nlev,islope REAL z4st,zdplanck(ngrid) @@ -106,5 +108,5 @@ REAL zkq(ngrid,nlay+1) REAL zcdv(ngrid),zcdh(ngrid) - REAL zcdv_true(ngrid),zcdh_true(ngrid) ! drag coeff are used by the LES to recompute u* and hfx + REAL, INTENT(IN) :: zcdv_true(ngrid),zcdh_true(ngrid) ! drag coeff are used by the LES to recompute u* and hfx REAL zu(ngrid,nlay),zv(ngrid,nlay) REAL zh(ngrid,nlay) @@ -136,7 +138,7 @@ REAL zdqsdif(ngrid) ! subtimestep pdqsdif for water ice REAL ztsrf(ngrid) ! temporary surface temperature in tsub - REAL zdtsrf(ngrid) ! surface temperature tendancy in tsub - REAL surf_h2o_lh(ngrid) ! Surface h2o latent heat flux - REAL zsurf_h2o_lh(ngrid) ! Tsub surface h2o latent heat flux + REAL zdtsrf(ngrid,nslope) ! surface temperature tendancy in tsub + REAL surf_h2o_lh(ngrid,nslope) ! Surface h2o latent heat flux + REAL zsurf_h2o_lh(ngrid,nslope) ! Tsub surface h2o latent heat flux c For latent heat release from ground water ice sublimation @@ -154,13 +156,15 @@ REAL qsat(ngrid) - REAL hdoflux(ngrid) ! value of vapour flux of HDO + REAL hdoflux(ngrid,nslope) ! value of vapour flux of HDO + REAL hdoflux_meshavg(ngrid) ! value of vapour flux of HDO ! REAL h2oflux(ngrid) ! value of vapour flux of H2O REAL old_h2o_vap(ngrid) ! traceur d'eau avant traitement + REAL saved_h2o_vap(ngrid) ! traceur d'eau avant traitement REAL kmixmin c Argument added for surface water ice budget: - REAL,INTENT(IN) :: watercap(ngrid) - REAL,INTENT(OUT) :: dwatercap_dif(ngrid) + REAL,INTENT(IN) :: watercap(ngrid,nslope) + REAL,INTENT(OUT) :: dwatercap_dif(ngrid,nslope) c Subtimestep to compute h2o latent heat flux: @@ -198,4 +202,16 @@ !!MARGAUX REAL DoH_vap(ngrid,nlay) +!! Sub-grid scale slopes + REAL :: pdqsdif_tmp(ngrid,nq) ! Temporary for dust lifting + REAL :: watercap_tmp(ngrid) + REAL :: zq_slope_vap(ngrid,nlay,nq,nslope) + REAL :: zq_tmp_vap(ngrid,nlay,nq) + REAL :: ptsrf_tmp(ngrid) + REAL :: pqsurf_tmp(ngrid,nq) + REAL :: pdqsdif_tmphdo(ngrid,nq) + REAL :: qsat_tmp(ngrid) + INTEGER :: indmax + + character*2 str2 c ** un petit test de coherence @@ -232,4 +248,8 @@ ENDIF + DO ig = 1,ngrid + ptsrf_tmp(ig) = ptsrf(ig,major_slope(ig)) + pqsurf_tmp(ig,:) = pqsurf(ig,:,major_slope(ig)) + ENDDO c----------------------------------------------------------------------- @@ -243,12 +263,13 @@ pdvdif(1:ngrid,1:nlay)=0 pdhdif(1:ngrid,1:nlay)=0 - pdtsrf(1:ngrid)=0 - zdtsrf(1:ngrid)=0 - surf_h2o_lh(1:ngrid)=0 - zsurf_h2o_lh(1:ngrid)=0 + pdtsrf(1:ngrid,1:nslope)=0 + zdtsrf(1:ngrid,1:nslope)=0 + surf_h2o_lh(1:ngrid,1:nslope)=0 + zsurf_h2o_lh(1:ngrid,1:nslope)=0 pdqdif(1:ngrid,1:nlay,1:nq)=0 - pdqsdif(1:ngrid,1:nq)=0 + pdqsdif(1:ngrid,1:nq,1:nslope)=0 + pdqsdif_tmp(1:ngrid,1:nq)=0 zdqsdif(1:ngrid)=0 - dwatercap_dif(1:ngrid)=0 + dwatercap_dif(1:ngrid,1:nslope)=0 c ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp @@ -275,5 +296,5 @@ ENDDO DO ig=1,ngrid - zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig)) + zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf_tmp(ig)) ENDDO @@ -354,4 +375,5 @@ ENDDO ENDDO + zq(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)+ & pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep @@ -366,6 +388,6 @@ c --------------------- - CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf,ph - & ,zcdv_true,zcdh_true) + CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf_tmp + & ,ph,zcdv_true,zcdh_true) zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1) @@ -383,5 +405,6 @@ DO ig=1,ngrid IF (zcdh_true(ig) .ne. 0.) THEN ! When Cd=Ch=0, u*=t*=0 - tstar(ig)=(ph(ig,1)-ptsrf(ig))*zcdh(ig)/ustar(ig) + tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig)) + & *zcdh(ig)/ustar(ig) ENDIF ENDDO @@ -391,5 +414,6 @@ zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)) - tstar(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:)) + tstar(:)=(ph(:,1)-ptsrf_tmp(:)) + & *zcdh_true(:)/sqrt(zcdv_true(:)) ENDIF @@ -454,7 +478,4 @@ ENDIF - - - c----------------------------------------------------------------------- c 4. inversion pour l'implicite sur u @@ -498,8 +519,4 @@ ENDDO - - - - c----------------------------------------------------------------------- c 5. inversion pour l'implicite sur v @@ -539,8 +556,4 @@ ENDDO ENDDO - - - - c----------------------------------------------------------------------- @@ -575,5 +588,7 @@ IF (tke_heat_flux .eq. 0.) THEN DO ig=1,ngrid - zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig) + indmax = major_slope(ig) + zdplanck(ig)=z4st*pemis(ig,indmax)*ptsrf(ig,indmax)* + & ptsrf(ig,indmax)*ptsrf(ig,indmax) ENDDO ELSE @@ -603,6 +618,8 @@ ENDIF + + DO ig=1,ngrid - + indmax = major_slope(ig) ! Initialization of z1 and zd, which do not depend on dmice @@ -649,7 +666,10 @@ c ---------------------------------------------------- - z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) - s +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep - z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig) + z1(ig)=pcapcal(ig,indmax)*ptsrf(ig,indmax) + s + cpp*zb(ig,1)*zc(ig,1) + s + zdplanck(ig)*ptsrf(ig,indmax) + s + pfluxsrf(ig,indmax)*ptimestep + z2(ig)= pcapcal(ig,indmax)+cpp*zb(ig,1)*(1.-zd(ig,1)) + s +zdplanck(ig) ztsrf2(ig)=z1(ig)/z2(ig) ! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop @@ -687,14 +707,35 @@ ENDDO ENDIF - + ENDDO!of Do j=1,XXX - + pdtsrf(ig,indmax)=(ztsrf2(ig)-ptsrf(ig,indmax))/ptimestep ENDDO !of Do ig=1,ngrid - - pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep DO ig=1,ngrid ! computing sensible heat flux (atm => surface) sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zhs(ig,1)-ztsrf2(ig)) ENDDO + +c Now implicit sheme on each sub-grid subslope: + IF (nslope.ne.1) then + DO islope=1,nslope + DO ig=1,ngrid + IF(islope.ne.major_slope(ig)) then + IF (tke_heat_flux .eq. 0.) THEN + zdplanck(ig)=z4st*pemis(ig,islope)*ptsrf(ig,islope)**3 + ELSE + zdplanck(ig) = 0. + ENDIF + z1(ig)=pcapcal(ig,islope)*ptsrf(ig,islope) + s + cpp*zb(ig,1)*zc(ig,1) + s + zdplanck(ig)*ptsrf(ig,islope) + s + pfluxsrf(ig,islope)*ptimestep + z2(ig)= pcapcal(ig,islope)+cpp*zb(ig,1)*(1.-zd(ig,1)) + s +zdplanck(ig) + ztsrf2(ig)=z1(ig)/z2(ig) + pdtsrf(ig,islope)=(ztsrf2(ig)-ptsrf(ig,islope))/ptimestep + ENDIF ! islope != indmax + ENDDO ! ig + ENDDO !islope + ENDIF !nslope.ne.1 c----------------------------------------------------------------------- @@ -727,9 +768,9 @@ do ig=1,ngrid c if(qsurf(ig,igcm_co2).lt.1) then - pdqsdif(ig,igcm_dust_mass) = + pdqsdif_tmp(ig,igcm_dust_mass) = & -alpha_lift(igcm_dust_mass) - pdqsdif(ig,igcm_dust_number) = + pdqsdif_tmp(ig,igcm_dust_number) = & -alpha_lift(igcm_dust_number) - pdqsdif(ig,igcm_dust_submicron) = + pdqsdif_tmp(ig,igcm_dust_submicron) = & -alpha_lift(igcm_dust_submicron) c end if @@ -739,8 +780,8 @@ !or 2 (injection in CL) do ig=1,ngrid - if(pqsurf(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2 - pdqsdif(ig,igcm_dust_mass) = + if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2 + pdqsdif_tmp(ig,igcm_dust_mass) = & -alpha_lift(igcm_dust_mass) - pdqsdif(ig,igcm_dust_number) = + pdqsdif_tmp(ig,igcm_dust_number) = & -alpha_lift(igcm_dust_number) end if @@ -748,33 +789,33 @@ elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface do ig=1,ngrid - if(pqsurf(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2 + if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2 IF((ti_injection_sol.LE.local_time(ig)).and. & (local_time(ig).LE.tf_injection_sol)) THEN if (rdstorm) then !Rocket dust storm scheme - pdqsdif(ig,igcm_stormdust_mass) = + pdqsdif_tmp(ig,igcm_stormdust_mass) = & -alpha_lift(igcm_stormdust_mass) & *dustliftday(ig) - pdqsdif(ig,igcm_stormdust_number) = + pdqsdif_tmp(ig,igcm_stormdust_number) = & -alpha_lift(igcm_stormdust_number) & *dustliftday(ig) - pdqsdif(ig,igcm_dust_mass)= 0. - pdqsdif(ig,igcm_dust_number)= 0. + pdqsdif_tmp(ig,igcm_dust_mass)= 0. + pdqsdif_tmp(ig,igcm_dust_number)= 0. else - pdqsdif(ig,igcm_dust_mass)= + pdqsdif_tmp(ig,igcm_dust_mass)= & -dustliftday(ig)* & alpha_lift(igcm_dust_mass) - pdqsdif(ig,igcm_dust_number)= + pdqsdif_tmp(ig,igcm_dust_number)= & -dustliftday(ig)* & alpha_lift(igcm_dust_number) endif if (submicron) then - pdqsdif(ig,igcm_dust_submicron) = 0. + pdqsdif_tmp(ig,igcm_dust_submicron) = 0. endif ! if (submicron) ELSE ! outside dust injection time frame - pdqsdif(ig,igcm_dust_mass)= 0. - pdqsdif(ig,igcm_dust_number)= 0. + pdqsdif_tmp(ig,igcm_dust_mass)= 0. + pdqsdif_tmp(ig,igcm_dust_number)= 0. if (rdstorm) then - pdqsdif(ig,igcm_stormdust_mass)= 0. - pdqsdif(ig,igcm_stormdust_number)= 0. + pdqsdif_tmp(ig,igcm_stormdust_mass)= 0. + pdqsdif_tmp(ig,igcm_stormdust_number)= 0. end if ENDIF @@ -785,5 +826,5 @@ else if (submicron) then do ig=1,ngrid - pdqsdif(ig,igcm_dust_submicron) = + pdqsdif_tmp(ig,igcm_dust_submicron) = & -alpha_lift(igcm_dust_submicron) end do @@ -791,10 +832,10 @@ #endif call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh, - & pqsurf(:,igcm_co2),pdqsdif) + & pqsurf_tmp(:,igcm_co2),pdqsdif_tmp) #ifndef MESOSCALE endif !doubleq.AND.submicron #endif else - pdqsdif(1:ngrid,1:nq) = 0. + pdqsdif_tmp(1:ngrid,1:nq) = 0. end if @@ -847,5 +888,5 @@ zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ $ zb(ig,2)*zc(ig,2) + - $ (-pdqsdif(ig,iq)) *ptimestep) *z1(ig) !tracer flux from surface + $ (-pdqsdif_tmp(ig,iq)) *ptimestep) *z1(ig) !tracer flux from surface ENDDO endif !((iq.eq.igcm_h2o_ice) @@ -857,4 +898,11 @@ ENDDO ENDDO + DO islope = 1,nslope + DO ig = 1,ngrid + pdqsdif(ig,iq,islope) = pdqsdif_tmp(ig,iq) + & * cos(pi*def_slope_mean(islope)/180.) + ENDDO + ENDDO + endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then enddo ! of do iq=1,nq @@ -867,11 +915,13 @@ c de decrire le flux de chaleur latente - do iq=1,nq if ((water).and.(iq.eq.igcm_h2o_vap)) then - + DO islope = 1,nslope DO ig=1,ngrid - zqsurf(ig)=pqsurf(ig,igcm_h2o_ice) + zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/ + & cos(pi*def_slope_mean(islope)/180.) + watercap_tmp(ig) = watercap(ig,islope)/ + & cos(pi*def_slope_mean(islope)/180.) ENDDO ! ig=1,ngrid @@ -879,5 +929,5 @@ c la subroutine est a la fin du fichier - call make_tsub(ngrid,pdtsrf,zqsurf, + call make_tsub(ngrid,pdtsrf(:,islope),zqsurf, & ptimestep,dtmax,watercaptag, & nsubtimestep) @@ -886,8 +936,10 @@ c initialization of ztsrf, which is surface temperature in c the subtimestep. + saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap) + DO ig=1,ngrid subtimestep = ptimestep/nsubtimestep(ig) - ztsrf(ig)=ptsrf(ig) ! +pdtsrf(ig)*subtimestep - + ztsrf(ig)=ptsrf(ig,islope) ! +pdtsrf(ig)*subtimestep + zq_tmp_vap(ig,:,:) =zq(ig,:,:) c Debut du sous pas de temps @@ -895,5 +947,4 @@ c C'est parti ! - zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) & /float(nsubtimestep(ig)) @@ -903,10 +954,10 @@ z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) - zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) + 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(ig,ilay,iq)+ + 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) @@ -914,9 +965,9 @@ z1(ig)=1./(za(ig,1)+zb(ig,1)+ $ zb(ig,2)*(1.-zd(ig,2))) - zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ + zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+ $ zb(ig,2)*zc(ig,2)) * z1(ig) call watersat(1,ztsrf(ig),pplev(ig,1),qsat(ig)) - old_h2o_vap(ig)=zq(ig,1,igcm_h2o_vap) + old_h2o_vap(ig)=zq_tmp_vap(ig,1,igcm_h2o_vap) zd(ig,1)=zb(ig,1)*z1(ig) zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig) @@ -925,4 +976,5 @@ & *(zq1temp(ig)-qsat(ig)) c write(*,*)'subliming more than available frost: qsurf!' + if(.not.watercaptag(ig)) then if ((-zdqsdif(ig)*subtimestep) @@ -935,5 +987,5 @@ c write(*,*)'flux vers le sol=',pdqsdif(ig,nq) z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2))) - zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+ + zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,igcm_h2o_vap)+ $ zb(ig,2)*zc(ig,2) + $ (-zdqsdif(ig)) *subtimestep) *z1(ig) @@ -944,5 +996,5 @@ c Starting upward calculations for water : c Actualisation de h2o_vap dans le premier niveau - zq(ig,1,igcm_h2o_vap)=zq1temp(ig) + zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig) c Take into account the H2O latent heat impact on the surface temperature @@ -950,22 +1002,23 @@ lh=(2834.3-0.28*(ztsrf(ig)-To)- & 0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3 - zdtsrf(ig)= zdqsdif(ig)*lh /pcapcal(ig) + zdtsrf(ig,islope)= zdqsdif(ig)*lh /pcapcal(ig,islope) endif ! (latentheat_surfwater) then DO ilay=2,nlay - zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) + zq_tmp_vap(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay) + & *zq_tmp_vap(ig,ilay-1,iq) ENDDO c Subtimestep water budget : - ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig) + zdtsrf(ig)) - & *subtimestep + ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope) + & + zdtsrf(ig,islope))*subtimestep zqsurf(ig)= zqsurf(ig)+( & zdqsdif(ig))*subtimestep c Monitoring instantaneous latent heat flux in W.m-2 : - zsurf_h2o_lh(ig) = zsurf_h2o_lh(ig)+ - & (zdtsrf(ig)*pcapcal(ig)) - & *subtimestep + zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+ + & (zdtsrf(ig,islope)*pcapcal(ig,islope)) + & *subtimestep c We ensure that surface temperature can't rise above the solid domain if there @@ -973,8 +1026,8 @@ if(zqsurf(ig) & +zdqsdif(ig)*subtimestep - & .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To - & zdtsrf(ig) = min(zdtsrf(ig),(To-ztsrf(ig))/subtimestep) ! ice melt case - - + & .gt.frost_albedo_threshold) then ! if there is still ice, T cannot exceed To + zdtsrf(ig,islope) = min(zdtsrf(ig,islope), + & (To-ztsrf(ig))/subtimestep) ! ice melt case + endif c Fin du sous pas de temps @@ -984,23 +1037,56 @@ c (btw could also compute the post timestep temp and ice c by simply adding the subtimestep trend instead of this) - surf_h2o_lh(ig)= zsurf_h2o_lh(ig)/ptimestep - pdtsrf(ig)= (ztsrf(ig) - - & ptsrf(ig))/ptimestep - pdqsdif(ig,igcm_h2o_ice)= - & (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice))/ptimestep + surf_h2o_lh(ig,islope)= zsurf_h2o_lh(ig,islope)/ptimestep + pdtsrf(ig,islope)= (ztsrf(ig) - + & ptsrf(ig,islope))/ptimestep + pdqsdif(ig,igcm_h2o_ice,islope)= + & (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice,islope)/ + & 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) if(watercaptag(ig)) then - if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep) - & .gt.(pqsurf(ig,igcm_h2o_ice))) then - dwatercap_dif(ig)=pdqsdif(ig,igcm_h2o_ice)+ - & pqsurf(ig,igcm_h2o_ice)/ptimestep - pdqsdif(ig,igcm_h2o_ice)= - & - pqsurf(ig,igcm_h2o_ice)/ptimestep + if ((-pdqsdif(ig,igcm_h2o_ice,islope)*ptimestep) + & .gt.(pqsurf(ig,igcm_h2o_ice,islope) + & /cos(pi*def_slope_mean(islope)/180.))) then + dwatercap_dif(ig,islope)= + & pdqsdif(ig,igcm_h2o_ice,islope)+ + & (pqsurf(ig,igcm_h2o_ice,islope) / + & cos(pi*def_slope_mean(islope)/180.))/ptimestep + pdqsdif(ig,igcm_h2o_ice,islope)= + & - (pqsurf(ig,igcm_h2o_ice,islope)/ + & cos(pi*def_slope_mean(islope)/180.))/ptimestep endif! ((-pdqsdif(ig)*ptimestep) endif !(watercaptag(ig)) then - + zq_slope_vap(ig,:,:,islope) = zq_tmp_vap(ig,:,:) ENDDO ! of DO ig=1,ngrid + ENDDO ! islope + +c Some grid box averages: interface with the atmosphere + DO ig = 1,ngrid + DO ilay = 1,nlay + zq(ig,ilay,iq) = 0. + DO islope = 1,nslope + zq(ig,ilay,iq) = zq(ig,ilay,iq) + + $ zq_slope_vap(ig,ilay,iq,islope) * + $ subslope_dist(ig,islope) + ENDDO + ENDDO + ENDDO + +! Recompute values in kg/m^2 slopped + DO ig = 1,ngrid + DO islope = 1,nslope + pdqsdif(ig,igcm_h2o_ice,islope) = + & pdqsdif(ig,igcm_h2o_ice,islope) + & * cos(pi*def_slope_mean(islope)/180.) + + dwatercap_dif(ig,islope) = + & dwatercap_dif(ig,islope) + & * cos(pi*def_slope_mean(islope)/180.) + ENDDO + ENDDO + END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap)) @@ -1031,9 +1117,29 @@ ENDDO ENDDO + hdoflux_meshavg(:) = 0. + DO islope = 1,nslope + + pdqsdif_tmphdo(:,:) = pdqsdif(:,:,islope) + & /cos(pi*def_slope_mean(islope)/180.) + + call watersat(ngrid,pdtsrf(:,islope)*ptimestep + + & ptsrf(:,islope),pplev(:,1),qsat_tmp) CALL hdo_surfex(ngrid,nlay,nq,ptimestep, - & zt,pplay,zq,pqsurf, - & old_h2o_vap,qsat,pdqsdif,dwatercap_dif, - & hdoflux) + & zt,pplay,zq,pqsurf(:,:,islope), + & saved_h2o_vap,qsat_tmp, + & pdqsdif_tmphdo, + & dwatercap_dif(:,islope)/cos(pi*def_slope_mean(islope)/180.), + & hdoflux(:,islope)) + + pdqsdif(:,:,islope) = pdqsdif_tmphdo(:,:) * + & cos(pi*def_slope_mean(islope)/180.) + DO ig = 1,ngrid + hdoflux_meshavg(ig) = hdoflux_meshavg(ig) + + & hdoflux(ig,islope)*subslope_dist(ig,islope) + + ENDDO !ig = 1,ngrid + ENDDO !islope = 1,nslope + DO ig=1,ngrid z1(ig)=1./(za(ig,1)+zb(ig,1)+ @@ -1041,5 +1147,5 @@ zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ $ zb(ig,2)*zc(ig,2) + - $ (-hdoflux(ig)) *ptimestep) *z1(ig) !tracer flux from surface + $ (-hdoflux_meshavg(ig)) *ptimestep) *z1(ig) !tracer flux from surface ENDDO @@ -1060,5 +1166,5 @@ call write_output("surf_h2o_lh", & "Ground ice latent heat flux", - & "W.m-2",surf_h2o_lh(:)) + & "W.m-2",surf_h2o_lh(:,iflat)) C Diagnostic output for HDO @@ -1066,5 +1172,5 @@ ! CALL write_output('hdoflux', ! & 'hdoflux', -! & ' ',hdoflux(:)) +! & ' ',hdoflux_meshavg(:)) ! CALL write_output('h2oflux', ! & 'h2oflux', @@ -1101,5 +1207,5 @@ WRITE(*,'(a10,3a15)') s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)' - PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1) + PRINT*,ptsrf(ngrid/2+1,:),ztsrf2(ngrid/2+1) DO ilev=1,nlay WRITE(*,'(4f15.7)')