Changeset 2896 for trunk/LMDZ.MARS


Ignore:
Timestamp:
Feb 14, 2023, 11:06:29 AM (22 months ago)
Author:
romain.vande
Message:

Mars PCM:
First modifications to introduce subslope parametrization in surface processes.
New variables are added in the new module comslope_mod to describe the subslopes:

_ nslope (number of subslope)
_ subslope_dist (distribution of the slopes)
_ def_slope (bound of the slope statistics / repartitions) etc..

These variables are added to the starfi file.
Other GCM variables, outputs and subroutines are not modified yet.
Only nslope=1 is accepted so far (ie. same as no subslope parametrization)

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/iostart.F90

    r2399 r2896  
    1414    INTEGER,SAVE :: idim6 ! "nlayer" dimension
    1515    INTEGER,SAVE :: idim7 ! "Time" dimension
     16    INTEGER,SAVE :: idim8 ! "nslope" dimension
     17    INTEGER,SAVE :: idim9 ! "inter slope" dimension
    1618    INTEGER,SAVE :: timeindex ! current time index (for time-dependent fields)
    1719    INTEGER,PARAMETER :: length=100 ! size of tab_cntrl array
    … …  
    466468  USE tracer_mod, only: nqmx
    467469  USE comsoil_h, only: nsoilmx
     470  USE comslope_mod, only: nslope
    468471  IMPLICIT NONE
    469472    CHARACTER(LEN=*),INTENT(IN) :: filename
    … …  
    554557      ENDIF
    555558
     559      ierr=NF90_DEF_DIM(nid_restart,"nslope",nslope,idim8)
     560      IF (ierr/=NF90_NOERR) THEN
     561        write(*,*)'phyredem: problem defining nslope dimension'
     562        write(*,*)trim(nf90_strerror(ierr))
     563        CALL ABORT
     564      ENDIF
     565
     566      ierr=NF90_DEF_DIM(nid_restart,"inter slope",nslope+1,idim9)
     567      IF (ierr/=NF90_NOERR) THEN
     568        write(*,*)'phyredem: problem defining inter slope dimension'
     569        write(*,*)trim(nf90_strerror(ierr))
     570        CALL ABORT
     571      ENDIF
     572
    556573      ierr=NF90_ENDDEF(nid_restart)
    557574      IF (ierr/=NF90_NOERR) THEN
    … …  
    632649  USE dimphy
    633650  USE comsoil_h, only: nsoilmx
     651  USE comslope_mod, ONLY: nslope
    634652  USE mod_grid_phy_lmdz
    635653  USE mod_phys_lmdz_para
    … …  
    819837        endif ! of if (.not.present(time))
    820838
     839      ELSE IF (field_size==nslope) THEN
     840        ! input is a 2D "subsurface field" array
     841        if (.not.present(time)) then ! for a time-independent field
     842          ierr = NF90_REDEF(nid_restart)
     843#ifdef NC_DOUBLE
     844          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
     845                            (/idim2,idim8/),nvarid)
     846#else
     847          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
     848                            (/idim2,idim8/),nvarid)
     849#endif
     850          if (ierr.ne.NF90_NOERR) then
     851            write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
     852            write(*,*)trim(nf90_strerror(ierr))
     853          endif
     854          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
     855          ierr = NF90_ENDDEF(nid_restart)
     856          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
     857        else
     858          ! check if the variable has already been defined:
     859          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
     860          if (ierr/=NF90_NOERR) then ! variable not found, define it
     861            ierr=NF90_REDEF(nid_restart)
     862#ifdef NC_DOUBLE
     863            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
     864                              (/idim2,idim8,idim7/),nvarid)
     865#else
     866            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
     867                              (/idim2,idim8,idim7/),nvarid)
     868#endif
     869           if (ierr.ne.NF90_NOERR) then
     870              write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
     871              write(*,*)trim(nf90_strerror(ierr))
     872            endif
     873            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
     874            ierr=NF90_ENDDEF(nid_restart)
     875          endif
     876          ! Write the variable
     877          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
     878                            start=(/1,1,timeindex/))
     879
     880        endif ! of if (.not.present(time))
     881
    821882      ELSE
    822883        PRINT *, "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name)
    … …  
    890951                    nf90_inq_dimid, nf90_inquire_dimension, NF90_INQ_VARID
    891952  USE comsoil_h, only: nsoilmx
     953  USE comslope_mod, only: nslope
    892954  USE mod_phys_lmdz_para, only: is_master
    893955  IMPLICIT NONE
    … …  
    9411003        ! We know it is an  "mlayer" kind of 1D array
    9421004        idim1d=idim3
     1005      ELSEIF (var_size==nslope+1) THEN
     1006        ! We know it is an  "inter slope" kind of 1D array
     1007        idim1d=idim9
    9431008      ELSE
    9441009        PRINT *, "put_var_rgen error : wrong dimension"

  • trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90

    r2892 r2896  
    1010subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
    1111                     day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf, &
    12                      tauscaling,totcloudfrac,wstar,watercap)
     12                     tauscaling,totcloudfrac,wstar,watercap,def_slope, &
     13                     def_slope_mean,subslope_dist)
    1314
    1415  use tracer_mod, only: noms ! tracer names
    … …  
    2526  USE ioipsl_getin_p_mod, ONLY : getin_p
    2627  use comsoil_h, only: flux_geo
     28  USE comslope_mod, ONLY: nslope, major_slope
    2729  implicit none
    2830 
    … …  
    6870  real,intent(out) :: wstar(ngrid) ! Max vertical velocity in thermals (m/s)
    6971  real,intent(out) :: watercap(ngrid) ! h2o_ice_cover
     72  real, intent(out) :: def_slope(nslope+1) !boundaries for bining of the slopes
     73  real, intent(out) :: def_slope_mean(nslope)
     74  real, intent(out) :: subslope_dist(ngrid,nslope) !undermesh statistics
    7075!======================================================================
    7176!  Local variables:
    … …  
    7479      real xmin,xmax ! to display min and max of a field
    7580!
    76       INTEGER ig,iq,lmax
     81      INTEGER ig,iq,lmax,islope
    7782      INTEGER nid, nvarid
    7883      INTEGER ierr, i, nsrf
    … …  
    106111      REAL :: watercaptag_tmp(ngrid)
    107112
     113! Sub-grid scale slopes
     114      LOGICAL :: startphy_slope !to be retrocompatible and add the nslope dimension
     115      REAL, ALLOCATABLE :: default_def_slope(:)
     116      REAL :: sum_dist
     117      REAL :: current_max !var to find max distrib slope
     118
    108119      CHARACTER(len=5) :: modname="phyetat0"
    109120
    … …  
    122133               p_omeg,p_g,p_mugaz,p_daysec,time0)
    123134endif ! of if (startphy_file)
     135
     136if(nslope.ne.1) then
     137  call abort_physic(modname, &
     138                "phyetat0: For now, nslope should be 1 (set in comslope_mod)",1)
     139endif
     140
     141allocate(default_def_slope(nslope+1))
     142!Sub-grid scale  subslopes
     143if (nslope.eq.7) then
     144  default_def_slope(1) = -43.
     145  default_def_slope(2) = -19.
     146  default_def_slope(3) = -9.
     147  default_def_slope(4) = -3.
     148  default_def_slope(5) = 3.
     149  default_def_slope(6) = 9.
     150  default_def_slope(7) = 19.
     151  default_def_slope(8) = 43.
     152elseif (nslope.eq.5) then
     153  default_def_slope(1) = -43.
     154  default_def_slope(2) = -9.
     155  default_def_slope(3) = -3.
     156  default_def_slope(4) = 3.
     157  default_def_slope(5) = 9.
     158  default_def_slope(6) = 43.
     159elseif (nslope.eq.1) then
     160  default_def_slope(1) = 0.
     161  default_def_slope(2) = 0.
     162endif
     163
     164if (startphy_file) then
     165 call get_var("def_slope",def_slope,found)
     166   if(.not.found) then
     167     startphy_slope=.false.
     168     write(*,*)'slope_settings: Problem while reading <def_slope>'
     169     write(*,*)'default def_slope will be used'
     170     do islope=1,nslope+1
     171       def_slope(islope) = default_def_slope(islope)
     172     enddo
     173     write(*,*)'computing corresponding distribution <subslope_dist>'
     174     write(*,*)'For now, woth nslope=1, subslope_dist is straigforward'
     175     write(*,*)'Later this operation will be done by newstart with a specific routine'
     176     subslope_dist(:,:)=1.
     177     !call subslope_mola(ngrid,nslope,def_slope,subslope_dist)
     178   else
     179     startphy_slope=.true.
     180     call get_field("subslope_dist",subslope_dist,found,indextime)
     181     if(.not.found) then
     182       write(*,*)'slope_settings: Problem while reading <subslope_dist>'
     183       write(*,*)'computing a new distribution'
     184       write(*,*)'For now, woth nslope=1, subslope_dist is straigforward'
     185       write(*,*)'Later this operation will be done by newstart with a specific routine'
     186       subslope_dist(:,:)=1.
     187     !call subslope_mola(ngrid,nslope,def_slope,subslope_dist)
     188     endif
     189    endif
     190else ! startphy_file 
     191 do islope=1,nslope+1
     192       def_slope(islope) = default_def_slope(islope)
     193     enddo
     194     write(*,*)'computing corresponding distribution <subslope_dist>'
     195     write(*,*)'For now, woth nslope=1, subslope_dist is straigforward'
     196     write(*,*)'Later this operation will be done by newstart with a specific routine'
     197     subslope_dist(:,:)=1.
     198     !call subslope_mola(ngrid,nslope,def_slope,subslope_dist)
     199endif   
     200
     201do islope=1,nslope
     202    def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
     203enddo
     204 
     205DO ig = 1,ngrid
     206  sum_dist = 0.
     207  DO islope = 1,nslope
     208     sum_dist = sum_dist + subslope_dist(ig,islope)
     209  ENDDO
     210  DO islope = 1,nslope
     211    subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
     212  ENDDO
     213ENDDO
     214
     215!Now determine the major subslope, ie. the maximal distribution
     216
     217DO ig=1,ngrid
     218  major_slope(ig)=1
     219  current_max=subslope_dist(ig,1)
     220  DO islope=2,nslope
     221    if(subslope_dist(ig,islope).GT.current_max) then
     222      major_slope(ig)=islope
     223      current_max=subslope_dist(ig,islope)
     224    ENDIF
     225  ENDDO
     226ENDDO
    124227
    125228if (startphy_file) then
    … …  
    637740write(*,*) "phyetat0: Surface water ice <watercap> range:", &
    638741            minval(watercap), maxval(watercap)
    639  
    640 
    641742
    642743if (startphy_file) then

  • trunk/LMDZ.MARS/libf/phymars/phyredem.F90

    r2887 r2896  
    77subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
    88                         phystep,day_ini,time,airefi, &
    9                          alb,ith)
     9                         alb,ith,def_slope, &
     10                         subslope_dist)
    1011! create physics restart file and write time-independent variables
    1112  use comsoil_h, only: inertiedat, volcapa, mlayer
    … …  
    2324  use comcstfi_h, only: g, mugaz, omeg, rad, rcp
    2425  use time_phylmdz_mod, only: daysec
     26  use comslope_mod, ONLY: nslope
    2527  implicit none
    2628 
    … …  
    3840  real,intent(in) :: alb(ngrid)
    3941  real,intent(in) :: ith(ngrid,nsoil)
     42  real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes
     43  real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics
    4044
    4145  real :: tab_cntrl(length) ! nb "length=100" defined in iostart module
    … …  
    144148
    145149  call put_field("watercaptag","Infinite water reservoir",watercaptag_tmp)
     150
     151  ! Sub grid scale slope parametrization
     152  call put_var("def_slope","slope criterium stages",def_slope)
     153  call put_field("subslope_dist","under mesh slope distribution",subslope_dist)
    146154
    147155  ! Close file

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r2887 r2896  
    104104      USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured
    105105      use ioipsl_getin_p_mod, only: getin_p
     106      use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
     107     &                        subslope_dist,iflat,sky_slope,
     108     &                        major_slope,compute_meshgridavg,
     109     &                        ini_comslope_h
    106110
    107111      IMPLICIT NONE
    … …  
    529533!$OMP THREADPRIVATE(check_physics_inputs,check_physics_outputs)     
    530534
     535c Sub-grid scale slopes
     536      integer :: islope
     537
    531538      logical :: write_restart
    532539
    … …  
    563570c        ~~~~~~~~~~~~
    564571#ifndef MESOSCALE
     572
     573         call ini_comslope_h(ngrid)
     574
    565575! GCM. Read netcdf initial physical parameters.
    566576         CALL phyetat0 ("startfi.nc",0,0,
    … …  
    569579     &         tsurf,tsoil,albedo,emis,
    570580     &         q2,qsurf,tauscaling,totcloudfrac,wstar,
    571      &         watercap)
     581     &         watercap,def_slope,def_slope_mean,subslope_dist)
     582
     583       DO islope=1,nslope
     584       sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2.
     585       END DO
     586
     587!     Sky view:
     588       DO islope=1,nslope
     589       sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2.
     590       END DO
     591!     Determine the 'flatest' slopes
     592       iflat = 1
     593       DO islope=2,nslope
     594         IF(abs(def_slope_mean(islope)).lt.
     595     &      abs(def_slope_mean(iflat)))THEN
     596           iflat = islope
     597         ENDIF
     598       ENDDO
     599       PRINT*,'Flat slope for islope = ',iflat
     600       PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
    572601
    573602#else
    … …  
    646675c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    647676         CALL surfini(ngrid,qsurf)
     677
    648678         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
    649679
    … …  
    714744     &                   nsoilmx,ngrid,nlayer,nq,
    715745     &                   ptimestep,pday,0.,cell_area,
    716      &                   albedodat,inertiedat)
     746     &                   albedodat,inertiedat,def_slope,
     747     &                   subslope_dist)
    717748          else
    718749             call physdem0("restartfi.nc",longitude,latitude,
    719750     &                   nsoilmx,ngrid,nlayer,nq,
    720751     &                   ptimestep,float(day_end),0.,cell_area,
    721      &                   albedodat,inertiedat)
     752     &                   albedodat,inertiedat,def_slope,
     753     &                   subslope_dist)
    722754          endif
    723755         endif
    … …  
    13611393         zcdh(:) = 0.
    13621394         zcdv(:) = 0.
     1395
    13631396         CALL vdifc(ngrid,nlayer,nq,zpopsk,
    13641397     $        ptimestep,capcal,lwrite,
    … …  
    13701403     &        zcondicea_co2microp,sensibFlux,
    13711404     &        dustliftday,local_time,watercap,dwatercap_dif)
     1405
    13721406          DO ig=1,ngrid
    13731407             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
    … …  
    18391873           call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
    18401874     &                zdqdev,zdqsdev)
    1841  
     1875
    18421876           if (dustbin.ge.1) then
    18431877              do iq=1,nq
    … …  
    18721906     &                pq,pdq,zdqsed,zdqssed,nq,
    18731907     &                tau,tauscaling)
     1908
    18741909c Flux at the surface of co2 ice computed in co2cloud microtimestep
    18751910           IF (rdstorm) THEN
    … …  
    19441979     $                   zdqcloud,zdqscloud,tau(:,1),qsurf(:,igcm_co2),
    19451980     $                   pu,pdu,pv,pdv,surfdust,surfice)
     1981
    19461982            endif ! of if (modulo(icount-1,ichemistry).eq.0)
    19471983
    … …  
    20402076     $              zdqssed_co2,zcondicea_co2microp,
    20412077     &              zdqsc)
     2078
    20422079         DO iq=1, nq
    20432080           DO ig=1,ngrid
Note: See TracChangeset for help on using the changeset viewer.