Index: trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90 (revision 2896) +++ trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90 (revision 2896) @@ -0,0 +1,88 @@ +MODULE comslope_mod +!======================================================================================================================! +! Subject: +!--------- +! Module used for the parametrization of subgrid scale slope +!----------------------------------------------------------------------------------------------------------------------! +! Reference: +!----------- +! Lange et al. (2023 in prep.), 'Introduction of Sub-Grid Scale Slope Microclimates in the Mars Planetary Climate Model', JGR +! +!======================================================================================================================! + +IMPLICIT NONE + +INTEGER, parameter :: nslope = 1 ! Number of slopes for the statistic (1) +INTEGER, SAVE :: iflat ! Number of slopes for the statistic (1) +REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: def_slope ! bound of the slope statistics / repartitions (°) +REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: def_slope_mean ! mean slope of each bin (°) +REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: sky_slope ! portion of the sky view by each slopes (1) +REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: subslope_dist ! distribution of the slopes (1) +INTEGER,SAVE,ALLOCATABLE,DIMENSION(:) :: major_slope ! Index of the subslope that occupies most of themesh (1) +!$OMP THREADPRIVATE(def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope) + +CONTAINS + SUBROUTINE ini_comslope_h(ngrid) + + IMPLICIT NONE + INTEGER, INTENT(IN) :: ngrid + + allocate(def_slope(nslope+1)) + allocate(def_slope_mean(nslope)) + allocate(sky_slope(nslope)) + allocate(subslope_dist(ngrid,nslope)) + allocate(major_slope(ngrid)) + END SUBROUTINE ini_comslope_h + + + SUBROUTINE end_comslope_h + + IMPLICIT NONE + + IF (allocated(def_slope)) deallocate(def_slope) + IF (allocated(def_slope_mean)) deallocate(def_slope_mean) + IF (allocated(sky_slope)) deallocate(sky_slope) + IF (allocated(subslope_dist)) deallocate(subslope_dist) + IF (allocated(major_slope)) deallocate(major_slope) + + END SUBROUTINE end_comslope_h + + SUBROUTINE compute_meshgridavg(ngrid,nq,albedo_slope,emis_slope,tsurf_slope,qsurf_slope,albedo_meshavg,emis_meshavg,tsurf_meshavg, qsurf_meshavg) + USE comcstfi_h, only: pi + IMPLICIT NONE + INTEGER, INTENT(IN) :: ngrid,nq ! # points on the physical grid, tracers (1) + REAL, INTENT(IN) :: albedo_slope(ngrid,2,nslope) ! albedo on each sub-grid slope (1) + REAL, INTENT(IN) :: emis_slope(ngrid,nslope) ! emissivity on each sub-grid slope (1) + REAL, INTENT(IN) :: tsurf_slope(ngrid,nslope) ! Surface Temperature on each sub-grid slope (K) + REAL, INTENT(IN) :: qsurf_slope(ngrid,nq,nslope) ! Surface Tracers on each sub-grid slope (kg/m^2 sloped) + REAL, INTENT(OUT) :: albedo_meshavg(ngrid,2) ! grid box average of the albedo (1) + REAL, INTENT(OUT) :: emis_meshavg(ngrid) ! grid box average of the emissivity (1) + REAL, INTENT(OUT) :: tsurf_meshavg(ngrid) ! grid box average of the surface temperature (K) + REAL, INTENT(OUT) :: qsurf_meshavg(ngrid,nq) ! grid box average of the surface tracers (kg/m^2 flat) +! Local + integer :: islope,ig,l,iq + +!------------------- + + albedo_meshavg(:,:) = 0. + emis_meshavg(:) = 0. + tsurf_meshavg(:) = 0. + qsurf_meshavg(:,:) = 0. + + DO ig = 1,ngrid + DO islope = 1,nslope + DO l = 1,2 + albedo_meshavg(ig,l) = albedo_meshavg(ig,l) + albedo_slope(ig,l,islope)*subslope_dist(ig,islope) + ENDDO + DO iq = 1,nq + qsurf_meshavg(ig,iq) = qsurf_meshavg(ig,iq) + qsurf_slope(ig,iq,islope)*subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.) + ENDDO + emis_meshavg(ig) = emis_meshavg(ig) + emis_slope(ig,islope)*subslope_dist(ig,islope) + tsurf_meshavg(ig) = tsurf_meshavg(ig) + (emis_slope(ig,islope)*tsurf_slope(ig,islope)**4)*subslope_dist(ig,islope) + ENDDO + tsurf_meshavg(ig) = (tsurf_meshavg(ig)/emis_meshavg(ig))**(0.25) + ENDDO + + END SUBROUTINE compute_meshgridavg + +END MODULE comslope_mod Index: trunk/LMDZ.MARS/libf/phymars/iostart.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/iostart.F90 (revision 2892) +++ trunk/LMDZ.MARS/libf/phymars/iostart.F90 (revision 2896) @@ -14,4 +14,6 @@ INTEGER,SAVE :: idim6 ! "nlayer" dimension INTEGER,SAVE :: idim7 ! "Time" dimension + INTEGER,SAVE :: idim8 ! "nslope" dimension + INTEGER,SAVE :: idim9 ! "inter slope" dimension INTEGER,SAVE :: timeindex ! current time index (for time-dependent fields) INTEGER,PARAMETER :: length=100 ! size of tab_cntrl array @@ -466,4 +468,5 @@ USE tracer_mod, only: nqmx USE comsoil_h, only: nsoilmx + USE comslope_mod, only: nslope IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: filename @@ -554,4 +557,18 @@ ENDIF + ierr=NF90_DEF_DIM(nid_restart,"nslope",nslope,idim8) + IF (ierr/=NF90_NOERR) THEN + write(*,*)'phyredem: problem defining nslope dimension' + write(*,*)trim(nf90_strerror(ierr)) + CALL ABORT + ENDIF + + ierr=NF90_DEF_DIM(nid_restart,"inter slope",nslope+1,idim9) + IF (ierr/=NF90_NOERR) THEN + write(*,*)'phyredem: problem defining inter slope dimension' + write(*,*)trim(nf90_strerror(ierr)) + CALL ABORT + ENDIF + ierr=NF90_ENDDEF(nid_restart) IF (ierr/=NF90_NOERR) THEN @@ -632,4 +649,5 @@ USE dimphy USE comsoil_h, only: nsoilmx + USE comslope_mod, ONLY: nslope USE mod_grid_phy_lmdz USE mod_phys_lmdz_para @@ -819,4 +837,47 @@ endif ! of if (.not.present(time)) + ELSE IF (field_size==nslope) THEN + ! input is a 2D "subsurface field" array + if (.not.present(time)) then ! for a time-independent field + ierr = NF90_REDEF(nid_restart) +#ifdef NC_DOUBLE + ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,& + (/idim2,idim8/),nvarid) +#else + ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,& + (/idim2,idim8/),nvarid) +#endif + if (ierr.ne.NF90_NOERR) then + write(*,*)"put_field_rgen error: failed to define"//trim(field_name) + write(*,*)trim(nf90_strerror(ierr)) + endif + IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title) + ierr = NF90_ENDDEF(nid_restart) + ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo) + else + ! check if the variable has already been defined: + ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid) + if (ierr/=NF90_NOERR) then ! variable not found, define it + ierr=NF90_REDEF(nid_restart) +#ifdef NC_DOUBLE + ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,& + (/idim2,idim8,idim7/),nvarid) +#else + ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,& + (/idim2,idim8,idim7/),nvarid) +#endif + if (ierr.ne.NF90_NOERR) then + write(*,*)"put_field_rgen error: failed to define"//trim(field_name) + write(*,*)trim(nf90_strerror(ierr)) + endif + IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title) + ierr=NF90_ENDDEF(nid_restart) + endif + ! Write the variable + ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,& + start=(/1,1,timeindex/)) + + endif ! of if (.not.present(time)) + ELSE PRINT *, "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name) @@ -890,4 +951,5 @@ nf90_inq_dimid, nf90_inquire_dimension, NF90_INQ_VARID USE comsoil_h, only: nsoilmx + USE comslope_mod, only: nslope USE mod_phys_lmdz_para, only: is_master IMPLICIT NONE @@ -941,4 +1003,7 @@ ! We know it is an "mlayer" kind of 1D array idim1d=idim3 + ELSEIF (var_size==nslope+1) THEN + ! We know it is an "inter slope" kind of 1D array + idim1d=idim9 ELSE PRINT *, "put_var_rgen error : wrong dimension" Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90 (revision 2892) +++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90 (revision 2896) @@ -10,5 +10,6 @@ subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, & day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf, & - tauscaling,totcloudfrac,wstar,watercap) + tauscaling,totcloudfrac,wstar,watercap,def_slope, & + def_slope_mean,subslope_dist) use tracer_mod, only: noms ! tracer names @@ -25,4 +26,5 @@ USE ioipsl_getin_p_mod, ONLY : getin_p use comsoil_h, only: flux_geo + USE comslope_mod, ONLY: nslope, major_slope implicit none @@ -68,4 +70,7 @@ real,intent(out) :: wstar(ngrid) ! Max vertical velocity in thermals (m/s) real,intent(out) :: watercap(ngrid) ! h2o_ice_cover + real, intent(out) :: def_slope(nslope+1) !boundaries for bining of the slopes + real, intent(out) :: def_slope_mean(nslope) + real, intent(out) :: subslope_dist(ngrid,nslope) !undermesh statistics !====================================================================== ! Local variables: @@ -74,5 +79,5 @@ real xmin,xmax ! to display min and max of a field ! - INTEGER ig,iq,lmax + INTEGER ig,iq,lmax,islope INTEGER nid, nvarid INTEGER ierr, i, nsrf @@ -106,4 +111,10 @@ REAL :: watercaptag_tmp(ngrid) +! Sub-grid scale slopes + LOGICAL :: startphy_slope !to be retrocompatible and add the nslope dimension + REAL, ALLOCATABLE :: default_def_slope(:) + REAL :: sum_dist + REAL :: current_max !var to find max distrib slope + CHARACTER(len=5) :: modname="phyetat0" @@ -122,4 +133,96 @@ p_omeg,p_g,p_mugaz,p_daysec,time0) endif ! of if (startphy_file) + +if(nslope.ne.1) then + call abort_physic(modname, & + "phyetat0: For now, nslope should be 1 (set in comslope_mod)",1) +endif + +allocate(default_def_slope(nslope+1)) +!Sub-grid scale subslopes +if (nslope.eq.7) then + default_def_slope(1) = -43. + default_def_slope(2) = -19. + default_def_slope(3) = -9. + default_def_slope(4) = -3. + default_def_slope(5) = 3. + default_def_slope(6) = 9. + default_def_slope(7) = 19. + default_def_slope(8) = 43. +elseif (nslope.eq.5) then + default_def_slope(1) = -43. + default_def_slope(2) = -9. + default_def_slope(3) = -3. + default_def_slope(4) = 3. + default_def_slope(5) = 9. + default_def_slope(6) = 43. +elseif (nslope.eq.1) then + default_def_slope(1) = 0. + default_def_slope(2) = 0. +endif + +if (startphy_file) then + call get_var("def_slope",def_slope,found) + if(.not.found) then + startphy_slope=.false. + write(*,*)'slope_settings: Problem while reading ' + write(*,*)'default def_slope will be used' + do islope=1,nslope+1 + def_slope(islope) = default_def_slope(islope) + enddo + write(*,*)'computing corresponding distribution ' + write(*,*)'For now, woth nslope=1, subslope_dist is straigforward' + write(*,*)'Later this operation will be done by newstart with a specific routine' + subslope_dist(:,:)=1. + !call subslope_mola(ngrid,nslope,def_slope,subslope_dist) + else + startphy_slope=.true. + call get_field("subslope_dist",subslope_dist,found,indextime) + if(.not.found) then + write(*,*)'slope_settings: Problem while reading ' + write(*,*)'computing a new distribution' + write(*,*)'For now, woth nslope=1, subslope_dist is straigforward' + write(*,*)'Later this operation will be done by newstart with a specific routine' + subslope_dist(:,:)=1. + !call subslope_mola(ngrid,nslope,def_slope,subslope_dist) + endif + endif +else ! startphy_file + do islope=1,nslope+1 + def_slope(islope) = default_def_slope(islope) + enddo + write(*,*)'computing corresponding distribution ' + write(*,*)'For now, woth nslope=1, subslope_dist is straigforward' + write(*,*)'Later this operation will be done by newstart with a specific routine' + subslope_dist(:,:)=1. + !call subslope_mola(ngrid,nslope,def_slope,subslope_dist) +endif + +do islope=1,nslope + def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2. +enddo + +DO ig = 1,ngrid + sum_dist = 0. + DO islope = 1,nslope + sum_dist = sum_dist + subslope_dist(ig,islope) + ENDDO + DO islope = 1,nslope + subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist + ENDDO +ENDDO + +!Now determine the major subslope, ie. the maximal distribution + +DO ig=1,ngrid + major_slope(ig)=1 + current_max=subslope_dist(ig,1) + DO islope=2,nslope + if(subslope_dist(ig,islope).GT.current_max) then + major_slope(ig)=islope + current_max=subslope_dist(ig,islope) + ENDIF + ENDDO +ENDDO if (startphy_file) then @@ -637,6 +740,4 @@ write(*,*) "phyetat0: Surface water ice range:", & minval(watercap), maxval(watercap) - - if (startphy_file) then Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/phyredem.F90 (revision 2892) +++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90 (revision 2896) @@ -7,5 +7,6 @@ subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, & phystep,day_ini,time,airefi, & - alb,ith) + alb,ith,def_slope, & + subslope_dist) ! create physics restart file and write time-independent variables use comsoil_h, only: inertiedat, volcapa, mlayer @@ -23,4 +24,5 @@ use comcstfi_h, only: g, mugaz, omeg, rad, rcp use time_phylmdz_mod, only: daysec + use comslope_mod, ONLY: nslope implicit none @@ -38,4 +40,6 @@ real,intent(in) :: alb(ngrid) real,intent(in) :: ith(ngrid,nsoil) + real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes + real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics real :: tab_cntrl(length) ! nb "length=100" defined in iostart module @@ -144,4 +148,8 @@ call put_field("watercaptag","Infinite water reservoir",watercaptag_tmp) + + ! Sub grid scale slope parametrization + call put_var("def_slope","slope criterium stages",def_slope) + call put_field("subslope_dist","under mesh slope distribution",subslope_dist) ! Close file Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 2892) +++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 2896) @@ -104,4 +104,8 @@ USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured use ioipsl_getin_p_mod, only: getin_p + use comslope_mod, ONLY: nslope,def_slope,def_slope_mean, + & subslope_dist,iflat,sky_slope, + & major_slope,compute_meshgridavg, + & ini_comslope_h IMPLICIT NONE @@ -529,4 +533,7 @@ !$OMP THREADPRIVATE(check_physics_inputs,check_physics_outputs) +c Sub-grid scale slopes + integer :: islope + logical :: write_restart @@ -563,4 +570,7 @@ c ~~~~~~~~~~~~ #ifndef MESOSCALE + + call ini_comslope_h(ngrid) + ! GCM. Read netcdf initial physical parameters. CALL phyetat0 ("startfi.nc",0,0, @@ -569,5 +579,24 @@ & tsurf,tsoil,albedo,emis, & q2,qsurf,tauscaling,totcloudfrac,wstar, - & watercap) + & watercap,def_slope,def_slope_mean,subslope_dist) + + DO islope=1,nslope + sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2. + END DO + +! Sky view: + DO islope=1,nslope + sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2. + END DO +! Determine the 'flatest' slopes + iflat = 1 + DO islope=2,nslope + IF(abs(def_slope_mean(islope)).lt. + & abs(def_slope_mean(iflat)))THEN + iflat = islope + ENDIF + ENDDO + PRINT*,'Flat slope for islope = ',iflat + PRINT*,'corresponding criterium = ',def_slope_mean(iflat) #else @@ -646,4 +675,5 @@ c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CALL surfini(ngrid,qsurf) + CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) @@ -714,10 +744,12 @@ & nsoilmx,ngrid,nlayer,nq, & ptimestep,pday,0.,cell_area, - & albedodat,inertiedat) + & albedodat,inertiedat,def_slope, + & subslope_dist) else call physdem0("restartfi.nc",longitude,latitude, & nsoilmx,ngrid,nlayer,nq, & ptimestep,float(day_end),0.,cell_area, - & albedodat,inertiedat) + & albedodat,inertiedat,def_slope, + & subslope_dist) endif endif @@ -1361,4 +1393,5 @@ zcdh(:) = 0. zcdv(:) = 0. + CALL vdifc(ngrid,nlayer,nq,zpopsk, $ ptimestep,capcal,lwrite, @@ -1370,4 +1403,5 @@ & zcondicea_co2microp,sensibFlux, & dustliftday,local_time,watercap,dwatercap_dif) + DO ig=1,ngrid zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig) @@ -1839,5 +1873,5 @@ call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2, & zdqdev,zdqsdev) - + if (dustbin.ge.1) then do iq=1,nq @@ -1872,4 +1906,5 @@ & pq,pdq,zdqsed,zdqssed,nq, & tau,tauscaling) + c Flux at the surface of co2 ice computed in co2cloud microtimestep IF (rdstorm) THEN @@ -1944,4 +1979,5 @@ $ zdqcloud,zdqscloud,tau(:,1),qsurf(:,igcm_co2), $ pu,pdu,pv,pdv,surfdust,surfice) + endif ! of if (modulo(icount-1,ichemistry).eq.0) @@ -2040,4 +2076,5 @@ $ zdqssed_co2,zcondicea_co2microp, & zdqsc) + DO iq=1, nq DO ig=1,ngrid