Ignore:
Timestamp:
Jan 25, 2024, 6:06:27 PM (11 months ago)
Author:
jbclement
Message:

Mars PCM:
Some changes to prepare the introduction of slopes in 1D: in particular, the subroutine "getslopes.F90" and "param_slope.F90" are now inside the module "slope_mod.F90" + Few small cleanings throughout the code.
JBC

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
2 deleted
7 edited

Legend:

Unmodified
Added
Removed

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

    r2909 r3183  
    33! Subject:
    44!---------
    5 !   Module used for the parametrization of subgrid scale slope 
     5!   Module used for the parametrization of subgrid scale slope
    66!----------------------------------------------------------------------------------------------------------------------!
    77! Reference:
    … …  
    1111!======================================================================================================================!
    1212
    13 IMPLICIT NONE
     13implicit none
    1414
    15 INTEGER, SAVE :: nslope                                ! Number of slopes for the statistic (1)
    16 INTEGER, SAVE :: iflat                                      ! Number of slopes for the statistic (1)
    17 REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: def_slope             ! bound of the slope statistics / repartitions (°)
    18 REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: def_slope_mean        ! mean slope of each bin (°)
    19 REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: sky_slope             ! portion of the sky view by each slopes (1)
    20 REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: subslope_dist       ! distribution of the slopes (1)
    21 INTEGER,SAVE,ALLOCATABLE,DIMENSION(:) :: major_slope        ! Index of the subslope that occupies most of themesh (1)
     15integer, save                              :: nslope         ! Number of slopes for the statistic (1)
     16integer, save                              :: iflat          ! Flat slope for islope (1)
     17real,    save, allocatable, dimension(:)   :: def_slope      ! Bound of the slope statistics / repartitions (°)
     18real,    save, allocatable, dimension(:)   :: def_slope_mean ! Mean slope of each bin (°)
     19real,    save, allocatable, dimension(:)   :: sky_slope      ! Portion of the sky view by each slopes (1)
     20real,    save, allocatable, dimension(:,:) :: subslope_dist  ! Distribution of the slopes (1)
     21integer, save, allocatable, dimension(:)   :: major_slope    ! Index of the subslope that occupies most of themesh (1)
    2222!$OMP THREADPRIVATE(nslope,def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope)
    2323
    24 CONTAINS
    25   SUBROUTINE ini_comslope_h(ngrid,nslope_in)
     24!=======================================================================
     25contains
     26!=======================================================================
    2627
    27   IMPLICIT NONE
    28   INTEGER, INTENT(IN) :: ngrid
    29   INTEGER, INTENT(IN) :: nslope_in
     28SUBROUTINE ini_comslope_h(ngrid,nslope_in)
    3029
    31     allocate(def_slope(nslope_in+1))
    32     allocate(def_slope_mean(nslope_in))
    33     allocate(sky_slope(nslope_in))
    34     allocate(subslope_dist(ngrid,nslope_in))
    35     allocate(major_slope(ngrid))
    36   END SUBROUTINE ini_comslope_h
     30implicit none
    3731
     32integer, intent(in) :: ngrid
     33integer, intent(in) :: nslope_in
    3834
    39   SUBROUTINE end_comslope_h
     35allocate(def_slope(nslope_in+1))
     36allocate(def_slope_mean(nslope_in))
     37allocate(sky_slope(nslope_in))
     38allocate(subslope_dist(ngrid,nslope_in))
     39allocate(major_slope(ngrid))
    4040
    41   IMPLICIT NONE
     41END SUBROUTINE ini_comslope_h
    4242
    43         IF (allocated(def_slope)) deallocate(def_slope)
    44         IF (allocated(def_slope_mean)) deallocate(def_slope_mean)
    45         IF (allocated(sky_slope)) deallocate(sky_slope)
    46         IF (allocated(subslope_dist)) deallocate(subslope_dist)
    47         IF (allocated(major_slope)) deallocate(major_slope)
     43!=======================================================================
    4844
    49   END SUBROUTINE end_comslope_h
     45SUBROUTINE end_comslope_h
    5046
    51   SUBROUTINE compute_meshgridavg(ngrid,nq,albedo_slope,emis_slope,tsurf_slope,qsurf_slope,albedo_meshavg,emis_meshavg,tsurf_meshavg, qsurf_meshavg)
    52   USE comcstfi_h, only:  pi
    53   IMPLICIT NONE
    54   INTEGER, INTENT(IN) :: ngrid,nq                  !  # points on the physical grid, tracers  (1)
    55   REAL, INTENT(IN) :: albedo_slope(ngrid,2,nslope) !  albedo on each sub-grid slope (1)
    56   REAL, INTENT(IN) :: emis_slope(ngrid,nslope)     !  emissivity on each sub-grid slope (1)
    57   REAL, INTENT(IN) :: tsurf_slope(ngrid,nslope)    !  Surface Temperature on each sub-grid slope (K)
    58   REAL, INTENT(IN) :: qsurf_slope(ngrid,nq,nslope) !  Surface Tracers on each sub-grid slope (kg/m^2 sloped)
    59   REAL, INTENT(OUT) :: albedo_meshavg(ngrid,2)     !  grid box average of the albedo (1)
    60   REAL, INTENT(OUT) :: emis_meshavg(ngrid)         !  grid box average of the emissivity (1)
    61   REAL, INTENT(OUT) :: tsurf_meshavg(ngrid)        !  grid box average of the surface temperature (K)
    62   REAL, INTENT(OUT) :: qsurf_meshavg(ngrid,nq)     !  grid box average of the surface tracers (kg/m^2 flat)
     47implicit none
     48
     49if (allocated(def_slope)) deallocate(def_slope)
     50if (allocated(def_slope_mean)) deallocate(def_slope_mean)
     51if (allocated(sky_slope)) deallocate(sky_slope)
     52if (allocated(subslope_dist)) deallocate(subslope_dist)
     53if (allocated(major_slope)) deallocate(major_slope)
     54
     55END SUBROUTINE end_comslope_h
     56
     57!=======================================================================
     58
     59SUBROUTINE compute_meshgridavg(ngrid,nq,albedo_slope,emis_slope,tsurf_slope,qsurf_slope,albedo_meshavg,emis_meshavg,tsurf_meshavg, qsurf_meshavg)
     60
     61use comcstfi_h, only: pi
     62
     63implicit none
     64
     65integer, intent(in) :: ngrid, nq                             ! # points on the physical grid, tracers (1)
     66real, dimension(ngrid,2,nslope),  intent(in) :: albedo_slope ! albedo on each sub-grid slope (1)
     67real, dimension(ngrid,nslope),    intent(in) :: emis_slope   ! emissivity on each sub-grid slope (1)
     68real, dimension(ngrid,nslope),    intent(in) :: tsurf_slope  ! Surface Temperature on each sub-grid slope (K)
     69real, dimension(ngrid,nq,nslope), intent(in) :: qsurf_slope  ! Surface Tracers on each sub-grid slope (kg/m^2 sloped)
     70real, dimension(ngrid,2),  intent(out) :: albedo_meshavg ! grid box average of the albedo (1)
     71real, dimension(ngrid),    intent(out) :: emis_meshavg   ! grid box average of the emissivity (1)
     72real, dimension(ngrid),    intent(out) :: tsurf_meshavg  ! grid box average of the surface temperature (K)
     73real, dimension(ngrid,nq), intent(out) :: qsurf_meshavg  ! grid box average of the surface tracers (kg/m^2 flat)
    6374! Local
    64   integer :: islope,ig,l,iq
     75integer :: islope, ig, l, iq
    6576
    6677!-------------------
    6778
    68       albedo_meshavg(:,:) = 0.
    69       emis_meshavg(:) = 0.
    70       tsurf_meshavg(:) = 0.
    71       qsurf_meshavg(:,:) = 0.
     79albedo_meshavg = 0.
     80emis_meshavg = 0.
     81tsurf_meshavg = 0.
     82qsurf_meshavg = 0.
    7283
    73   if(nslope.eq.1) then
    74       albedo_meshavg(:,:) = albedo_slope(:,:,1)
    75       emis_meshavg(:) = emis_slope(:,1)
    76       tsurf_meshavg(:) = tsurf_slope(:,1)
    77       qsurf_meshavg(:,:) = qsurf_slope(:,:,1)
    78   else
     84if (nslope == 1) then
     85    albedo_meshavg = albedo_slope(:,:,1)
     86    emis_meshavg = emis_slope(:,1)
     87    tsurf_meshavg = tsurf_slope(:,1)
     88    qsurf_meshavg = qsurf_slope(:,:,1)
     89else
     90    do ig = 1,ngrid
     91        do islope = 1,nslope
     92            do l = 1,2
     93                albedo_meshavg(ig,l) = albedo_meshavg(ig,l) + albedo_slope(ig,l,islope)*subslope_dist(ig,islope)
     94            enddo
     95            do iq = 1,nq
     96                qsurf_meshavg(ig,iq) = qsurf_meshavg(ig,iq) + qsurf_slope(ig,iq,islope)*subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
     97            enddo
     98            emis_meshavg(ig) = emis_meshavg(ig) + emis_slope(ig,islope)*subslope_dist(ig,islope)
     99            tsurf_meshavg(ig) = tsurf_meshavg(ig) + (emis_slope(ig,islope)*tsurf_slope(ig,islope)**4)*subslope_dist(ig,islope)
     100        enddo
     101        tsurf_meshavg(ig) = (tsurf_meshavg(ig)/emis_meshavg(ig))**(0.25)
     102    enddo
     103endif
    79104
    80   DO ig = 1,ngrid
    81       DO islope = 1,nslope
    82           DO l = 1,2
    83             albedo_meshavg(ig,l) = albedo_meshavg(ig,l) + albedo_slope(ig,l,islope)*subslope_dist(ig,islope)
    84           ENDDO
    85           DO iq = 1,nq
    86              qsurf_meshavg(ig,iq) = qsurf_meshavg(ig,iq) + qsurf_slope(ig,iq,islope)*subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
    87           ENDDO
    88           emis_meshavg(ig) = emis_meshavg(ig) + emis_slope(ig,islope)*subslope_dist(ig,islope)
    89           tsurf_meshavg(ig) = tsurf_meshavg(ig) + (emis_slope(ig,islope)*tsurf_slope(ig,islope)**4)*subslope_dist(ig,islope)
    90       ENDDO
    91       tsurf_meshavg(ig) = (tsurf_meshavg(ig)/emis_meshavg(ig))**(0.25)
    92   ENDDO
    93 
    94   endif
    95 
    96   END SUBROUTINE compute_meshgridavg
     105END SUBROUTINE compute_meshgridavg
    97106
    98107END MODULE comslope_mod

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

    r3126 r3183  
    1 module comsoil_h
     1MODULE comsoil_h
    22
    33implicit none
     4
    45! nsoilmx : number of subterranean layers
    5   integer, parameter :: nsoilmx = 57
     6integer, parameter :: nsoilmx = 57
    67
    7   real,save,allocatable,dimension(:) :: layer      ! soil layer depths
    8   real,save,allocatable,dimension(:) :: mlayer     ! soil mid-layer depths
    9   real,save,allocatable,dimension(:,:) :: inertiedat ! soil thermal inertia for present climate
    10   real,save,allocatable,dimension(:,:,:) :: inertiesoil ! soil thermal inertia
    11   real,save :: volcapa    ! soil volumetric heat capacity
    12        ! NB: volcapa is read fromn control(35) from physicq start file
    13        !     in physdem (or set via tabfi, or initialized in
    14        !                 soil_settings.F)
     8real, save, allocatable, dimension(:)     :: layer       ! soil layer depths
     9real, save, allocatable, dimension(:)     :: mlayer      ! soil mid-layer depths
     10real, save, allocatable, dimension(:,:)   :: inertiedat  ! soil thermal inertia for present climate
     11real, save, allocatable, dimension(:,:,:) :: inertiesoil ! soil thermal inertia
     12real, save                                :: volcapa     ! soil volumetric heat capacity
     13! NB: volcapa is read from control(35) from physics start file
     14!     in physdem (or set via tabfi, or initialized in soil_settings.F)
    1515
    1616!$OMP THREADPRIVATE(layer,mlayer,inertiedat,inertiesoil,volcapa)
    1717
    18   ! variables (FC: built in firstcall in soil.F)
    19   REAL,SAVE,ALLOCATABLE :: tsoil(:,:,:)       ! sub-surface temperatures (K)
    20   real,save,allocatable :: mthermdiff(:,:,:) ! (FC) mid-layer thermal diffusivity
    21   real,save,allocatable :: thermdiff(:,:,:)   ! (FC) inter-layer thermal diffusivity
    22   real,save,allocatable :: coefq(:)           ! (FC) q_{k+1/2} coefficients
    23   real,save,allocatable :: coefd(:,:,:)       ! (FC) d_k coefficients
    24   real,save,allocatable :: alph(:,:,:)        ! (FC) alpha_k coefficients
    25   real,save,allocatable :: beta(:,:,:)        ! beta_k coefficients
    26   real,save :: mu
    27   real,save,allocatable :: flux_geo(:,:)      ! Geothermal Flux (W/m^2)
     18! Variables (FC: built in firstcall in soil.F)
     19real, save, allocatable, dimension(:,:,:) :: tsoil      ! sub-surface temperatures (K)
     20real, save, allocatable, dimension(:,:,:) :: mthermdiff ! (FC) mid-layer thermal diffusivity
     21real, save, allocatable, dimension(:,:,:) :: thermdiff  ! (FC) inter-layer thermal diffusivity
     22real, save, allocatable, dimension(:)     :: coefq      ! (FC) q_{k+1/2} coefficients
     23real, save, allocatable, dimension(:,:,:) :: coefd      ! (FC) d_k coefficients
     24real, save, allocatable, dimension(:,:,:) :: alph       ! (FC) alpha_k coefficients
     25real, save, allocatable, dimension(:,:,:) :: beta       ! beta_k coefficients
     26real, save, allocatable, dimension(:,:)   :: flux_geo   ! Geothermal Flux (W/m^2)
     27real, save :: mu
    2828
    2929!$OMP THREADPRIVATE(tsoil,mthermdiff,thermdiff,coefq,coefd,alph,beta,mu,flux_geo)
    3030
     31! Subsurface tracers:
     32logical, save                                  :: adsorption_soil ! boolean to call adosrption (or not)
     33real,    save                                  :: choice_ads      ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
     34real,    save, allocatable, dimension(:,:,:,:) :: qsoil           ! subsurface tracers (kg/m^3 of regol)
     35integer, parameter :: nqsoil = 3                                  ! number of subsurface tracers, only three when working with water
     36integer, parameter :: igcm_h2o_vap_soil = 1
     37integer, parameter :: igcm_h2o_ice_soil = 2
     38integer, parameter :: igcm_h2o_vap_ads  = 3
    3139
    32 ! Subsurface tracers:
    33   logical,save :: adsorption_soil             ! boolean to call adosrption (or not)
    34   real,save :: choice_ads                     ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
    35   integer, parameter :: nqsoil = 3            ! number of subsurface tracers, only three when working with water
    36   real,save,allocatable :: qsoil(:,:,:,:)     ! subsurface tracers (kg/m^3 of regol)
    37   integer, parameter :: igcm_h2o_vap_soil = 1
    38   integer, parameter :: igcm_h2o_ice_soil = 2
    39   integer, parameter :: igcm_h2o_vap_ads  = 3
    4040!$OMP THREADPRIVATE(adsorption_soil,qsoil,choice_ads)
    4141
     42!=======================================================================
     43contains
     44!=======================================================================
    4245
    43 contains
     46subroutine ini_comsoil_h(ngrid,nslope)
    4447
    45   subroutine ini_comsoil_h(ngrid,nslope)
    46  
    47   implicit none
    48   integer,intent(in) :: ngrid ! number of atmospheric columns
    49   integer,intent(in) :: nslope ! number of sub grid slopes
    50     allocate(layer(nsoilmx)) !soil layer depths
    51     allocate(mlayer(0:nsoilmx-1)) ! soil mid-layer depths
    52     allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia for present climate
    53     allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
    54     allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
    55     allocate(mthermdiff(ngrid,0:nsoilmx-1,nslope))
    56     allocate(thermdiff(ngrid,nsoilmx-1,nslope))
    57     allocate(coefq(0:nsoilmx-1))
    58     allocate(coefd(ngrid,nsoilmx-1,nslope))
    59     allocate(alph(ngrid,nsoilmx-1,nslope))
    60     allocate(beta(ngrid,nsoilmx-1,nslope))
    61     allocate(flux_geo(ngrid,nslope))
    62     allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
    63  
    64   end subroutine ini_comsoil_h
     48implicit none
    6549
     50integer,intent(in) :: ngrid  ! number of atmospheric columns
     51integer,intent(in) :: nslope ! number of sub grid slopes
    6652
    67   subroutine end_comsoil_h
     53allocate(layer(nsoilmx)) !soil layer depths
     54allocate(mlayer(0:nsoilmx - 1)) ! soil mid-layer depths
     55allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia for present climate
     56allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
     57allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
     58allocate(mthermdiff(ngrid,0:nsoilmx - 1,nslope))
     59allocate(thermdiff(ngrid,nsoilmx - 1,nslope))
     60allocate(coefq(0:nsoilmx - 1))
     61allocate(coefd(ngrid,nsoilmx - 1,nslope))
     62allocate(alph(ngrid,nsoilmx - 1,nslope))
     63allocate(beta(ngrid,nsoilmx - 1,nslope))
     64allocate(flux_geo(ngrid,nslope))
     65allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
    6866
    69   implicit none
     67END SUBROUTINE ini_comsoil_h
    7068
    71     if (allocated(layer)) deallocate(layer)
    72     if (allocated(mlayer)) deallocate(mlayer)
    73     if (allocated(inertiedat)) deallocate(inertiedat)
    74     if (allocated(inertiesoil)) deallocate(inertiesoil)
    75     if (allocated(tsoil)) deallocate(tsoil)
    76     if (allocated(mthermdiff)) deallocate(mthermdiff)
    77     if (allocated(thermdiff)) deallocate(thermdiff)
    78     if (allocated(coefq)) deallocate(coefq)
    79     if (allocated(coefd)) deallocate(coefd)
    80     if (allocated(alph)) deallocate(alph)
    81     if (allocated(beta)) deallocate(beta)
    82     if (allocated(flux_geo)) deallocate(flux_geo)
    83     if (allocated(qsoil))  deallocate(qsoil)
    84   end subroutine end_comsoil_h
     69!=======================================================================
     70SUBROUTINE end_comsoil_h
    8571
    86   subroutine ini_comsoil_h_slope_var(ngrid,nslope)
    87  
    88   implicit none
    89   integer,intent(in) :: ngrid ! number of atmospheric columns
    90   integer,intent(in) :: nslope ! number of sub grid slopes
    91  
    92     allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
    93     allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
    94     allocate(mthermdiff(ngrid,0:nsoilmx-1,nslope))
    95     allocate(thermdiff(ngrid,nsoilmx-1,nslope))
    96     allocate(coefd(ngrid,nsoilmx-1,nslope))
    97     allocate(alph(ngrid,nsoilmx-1,nslope))
    98     allocate(beta(ngrid,nsoilmx-1,nslope))
    99     allocate(flux_geo(ngrid,nslope))
    100     allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
    101  
    102   end subroutine ini_comsoil_h_slope_var
     72implicit none
    10373
     74if (allocated(layer)) deallocate(layer)
     75if (allocated(mlayer)) deallocate(mlayer)
     76if (allocated(inertiedat)) deallocate(inertiedat)
     77if (allocated(inertiesoil)) deallocate(inertiesoil)
     78if (allocated(tsoil)) deallocate(tsoil)
     79if (allocated(mthermdiff)) deallocate(mthermdiff)
     80if (allocated(thermdiff)) deallocate(thermdiff)
     81if (allocated(coefq)) deallocate(coefq)
     82if (allocated(coefd)) deallocate(coefd)
     83if (allocated(alph)) deallocate(alph)
     84if (allocated(beta)) deallocate(beta)
     85if (allocated(flux_geo)) deallocate(flux_geo)
     86if (allocated(qsoil))  deallocate(qsoil)
    10487
    105   subroutine end_comsoil_h_slope_var
     88END SUBROUTINE end_comsoil_h
    10689
    107   implicit none
     90!=======================================================================
     91SUBROUTINE ini_comsoil_h_slope_var(ngrid,nslope)
    10892
    109     if (allocated(tsoil)) deallocate(tsoil)
    110     if (allocated(inertiesoil)) deallocate(inertiesoil)
    111     if (allocated(mthermdiff)) deallocate(mthermdiff)
    112     if (allocated(thermdiff)) deallocate(thermdiff)
    113     if (allocated(coefd)) deallocate(coefd)
    114     if (allocated(alph)) deallocate(alph)
    115     if (allocated(beta)) deallocate(beta)
    116     if (allocated(flux_geo)) deallocate(flux_geo)
    117     if (allocated(qsoil))  deallocate(qsoil)
     93implicit none
    11894
    119   end subroutine end_comsoil_h_slope_var
     95integer,intent(in) :: ngrid  ! number of atmospheric columns
     96integer,intent(in) :: nslope ! number of sub grid slopes
    12097
    121 end module comsoil_h
     98allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
     99allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
     100allocate(mthermdiff(ngrid,0:nsoilmx - 1,nslope))
     101allocate(thermdiff(ngrid,nsoilmx - 1,nslope))
     102allocate(coefd(ngrid,nsoilmx - 1,nslope))
     103allocate(alph(ngrid,nsoilmx - 1,nslope))
     104allocate(beta(ngrid,nsoilmx - 1,nslope))
     105allocate(flux_geo(ngrid,nslope))
     106allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
     107
     108END SUBROUTINE ini_comsoil_h_slope_var
     109
     110!=======================================================================
     111SUBROUTINE end_comsoil_h_slope_var
     112
     113implicit none
     114
     115if (allocated(tsoil)) deallocate(tsoil)
     116if (allocated(inertiesoil)) deallocate(inertiesoil)
     117if (allocated(mthermdiff)) deallocate(mthermdiff)
     118if (allocated(thermdiff)) deallocate(thermdiff)
     119if (allocated(coefd)) deallocate(coefd)
     120if (allocated(alph)) deallocate(alph)
     121if (allocated(beta)) deallocate(beta)
     122if (allocated(flux_geo)) deallocate(flux_geo)
     123if (allocated(qsoil))  deallocate(qsoil)
     124
     125END SUBROUTINE end_comsoil_h_slope_var
     126
     127END MODULE comsoil_h

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

    r3179 r3183  
    4141use conf_phys_mod,            only: conf_phys
    4242use paleoclimate_mod,         only: paleoclimate
     43use comslope_mod,             only: nslope, subslope_dist, ini_comslope_h, end_comslope_h
    4344! Mostly for XIOS outputs:
    4445use mod_const_mpi,            only: COMM_LMDZ
    … …  
    208209allocate(nqfils(nqtot))
    209210nqperes = 0
    210 nqfils(:) = 0
     211nqfils = 0
    211212do iq = 1,nqtot
    212213    if (tnom_transp(iq) == 'air') then
    … …  
    373374longitude = longitude*pi/180.
    374375
    375 ! Some initializations (some of which have already been
    376 ! done above!) and loads parameters set in callphys.def
    377 ! and allocates some arrays
     376! Some initializations (some of which have already been done above!)
     377! and loads parameters set in callphys.def and allocates some arrays
    378378! Mars possible matter with dttestphys in input and include!!!
    379379! Initializations below should mimick what is done in iniphysiq for 3D GCM
    … …  
    383383call init_geometry_cell_area_for_outputs(1,cell_area)
    384384! Ehouarn: init_vertial_layers called later (because disvert not called yet)
    385 !      call init_vertical_layers(nlayer,preff,scaleheight,
    386 !     &                      ap,bp,aps,bps,presnivs,pseudoalt)
     385! call init_vertical_layers(nlayer,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
    387386call init_dimphy(1,nlayer) ! Initialize dimphy module
    388 call phys_state_var_init(1,llm,nq,tname,day0,dayn,time,daysec,dttestphys,rad,g,r,cpp,nqperes,nqfils)! MVals: variables isotopes
     387call phys_state_var_init(1,llm,nq,tname,day0,dayn,time,daysec,dttestphys,rad,g,r,cpp,nqperes,nqfils) ! MVals: variables isotopes
    389388call ini_fillgeom(1,latitude,longitude,(/1.0/))
    390389call conf_phys(1,llm,nq)
    … …  
    461460
    462461! For the non-orographic gravity waves scheme
    463 du_nonoro_gwd(1,:)=0
    464 dv_nonoro_gwd(1,:)=0
     462du_nonoro_gwd(1,:) = 0
     463dv_nonoro_gwd(1,:) = 0
    465464
    466465! For the slope wind scheme
    … …  
    559558if (.not. therestart1D) then
    560559    tsurf(:,1) = tmp2(0)
    561     temp(:) = tmp2(1:)
     560    temp = tmp2(1:)
    562561else
    563562    read(3,*) header, (tsurf(1,:), j = 1,size(tsurf,2)), (temp(ilayer), ilayer = 1,nlayer)

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90

    r3179 r3183  
    255255    ! Compute winds and temperature for next time step
    256256    ! ------------------------------------------------
    257     u(:) = u(:) + dttestphys*du(:)
    258     v(:) = v(:) + dttestphys*dv(:)
    259     temp(:) = temp(:) + dttestphys*dtemp(:)
     257    u = u + dttestphys*du
     258    v = v + dttestphys*dv
     259    temp = temp + dttestphys*dtemp
    260260
    261261    ! Compute pressure for next time step

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

    r3174 r3183  
    146146
    147147if (startphy_file) then
    148  call get_var("def_slope",def_slope,found)
    149    if(.not.found) then
    150      startphy_slope=.false.
    151      write(*,*)'slope_settings: Problem while reading <def_slope>'
    152      write(*,*)'Checking that nslope=1'
    153      if(nslope.eq.1) then
    154       write(*,*)'We take default def_slope and subslope dist'
    155       allocate(default_def_slope(nslope+1))
    156       default_def_slope(1) = -90.
    157       default_def_slope(2) = 90.
    158       subslope_dist(:,:)=1.
    159       def_slope(:)=default_def_slope(:)
    160      else
    161        write(*,*)'Variable def_slope is not present in the start and'
    162        write(*,*)'you are trying to run with nslope!=1'
    163        write(*,*)'This is not possible, you have to go through newstart'
    164      endif
     148   call get_var("def_slope",def_slope,found)
     149   if (.not. found) then
     150       startphy_slope=.false.
     151       write(*,*)'slope_settings: Problem while reading <def_slope>'
     152       write(*,*)'Checking that nslope=1'
     153       if (nslope == 1) then
     154           write(*,*)'We take default def_slope and subslope dist'
     155           allocate(default_def_slope(nslope + 1))
     156           default_def_slope(1) = -90.
     157           default_def_slope(2) = 90.
     158           subslope_dist = 1.
     159           def_slope = default_def_slope
     160       else
     161           write(*,*)'Variable def_slope is not present in the start and'
     162           write(*,*)'you are trying to run with nslope!=1'
     163           write(*,*)'This is not possible, you have to go through newstart'
     164       endif
    165165   else
    166      startphy_slope=.true.
    167      call get_field("subslope_dist",subslope_dist,found,indextime)
    168      if(.not.found) then
    169        write(*,*)'slope_settings: Problem while reading <subslope_dist>'
    170        write(*,*)'We have to abort.'
    171        write(*,*)'Please check that nslope is coherent, you can modify this parameter with newstart'
    172        call abort_physic(modname, &
    173                 "phyetat0: Failed loading <subslope_dist>",1)
    174      endif
     166       startphy_slope=.true.
     167       call get_field("subslope_dist",subslope_dist,found,indextime)
     168       if (.not. found) then
     169           write(*,*)'slope_settings: Problem while reading <subslope_dist>'
     170           write(*,*)'We have to abort.'
     171           write(*,*)'Please check that nslope is coherent, you can modify this parameter with newstart'
     172           call abort_physic(modname,"phyetat0: Failed loading <subslope_dist>",1)
     173       endif
     174   endif
     175else ! (no startphy_file)
     176    if (nslope == 1) then
     177        allocate(default_def_slope(2))
     178        default_def_slope = 0.
     179        subslope_dist = 1.
     180        def_slope = default_def_slope
     181    else
     182        write(*,*)'Without startfi, lets first run with nslope=1'
     183        call abort_physic(modname,"phyetat0: No startfi and nslope!=1",1)
    175184    endif
    176 else ! (no startphy_file)
    177      if(nslope.eq.1) then
    178       allocate(default_def_slope(nslope+1))
    179       default_def_slope(1) = 0.
    180       default_def_slope(2) = 0.
    181       subslope_dist(:,:)=1.
    182       def_slope(:)=default_def_slope(:)
    183     else
    184      write(*,*)'Without starfi, lets first run with nslope=1'
    185        call abort_physic(modname, &
    186                 "phyetat0: No startfi and nslope!=1",1)
    187    endif
    188185endif   
    189186
    190 do islope=1,nslope
    191     def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
     187do islope = 1,nslope
     188    def_slope_mean(islope) = (def_slope(islope) + def_slope(islope + 1))/2.
    192189enddo
    193190 
    194191DO ig = 1,ngrid
    195   sum_dist = 0.
    196   DO islope = 1,nslope
    197      sum_dist = sum_dist + subslope_dist(ig,islope)
    198   ENDDO
    199   DO islope = 1,nslope
    200     subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
    201   ENDDO
     192    sum_dist = 0.
     193    DO islope = 1,nslope
     194        sum_dist = sum_dist + subslope_dist(ig,islope)
     195    ENDDO
     196    DO islope = 1,nslope
     197        subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
     198    ENDDO
    202199ENDDO
    203200

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

    r3167 r3183  
    6565      use nirdata_mod, only: NIR_leedat
    6666      use nirco2abs_mod, only: nirco2abs
    67       use slope_mod, only: theta_sl, psi_sl
     67      use slope_mod, only: theta_sl, psi_sl, getslopes, param_slope
    6868      use conc_mod, only: rnew, cpnew, mmean
    6969      use time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec

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

    r2578 r3183  
    1 module slope_mod
    2 
    3 implicit none
    4 
    5   real,save,allocatable :: theta_sl(:) ! slope angle versus horizontal (deg)
    6   real,save,allocatable :: psi_sl(:)   ! slope orientation (deg)
     1MODULE slope_mod
     2
     3implicit none
     4
     5real, save, dimension(:), allocatable :: theta_sl ! slope angle versus horizontal (deg)
     6real, save, dimension(:), allocatable :: psi_sl   ! slope orientation (deg)
    77
    88!$OMP THREADPRIVATE(theta_sl,psi_sl)
    99
     10!=======================================================================
    1011contains
    11 
    12   subroutine ini_slope_mod(ngrid)
    13  
    14   implicit none
    15   integer,intent(in) :: ngrid ! number of atmospheric columns
    16  
    17   allocate(theta_sl(ngrid))
    18   allocate(psi_sl(ngrid))
    19  
    20   end subroutine ini_slope_mod
    21 
    22 
    23   subroutine end_slope_mod
    24 
    25   implicit none
    26 
    27   if (allocated(theta_sl)) deallocate(theta_sl)
    28   if (allocated(psi_sl)) deallocate(psi_sl)
    29 
    30   end subroutine end_slope_mod
    31  
    32 end module slope_mod
     12!=======================================================================
     13
     14SUBROUTINE getslopes(ngrid,geopot)
     15
     16use geometry_mod,       only: longitude, latitude ! in radians
     17use comcstfi_h,         only: g, rad, pi
     18use mod_phys_lmdz_para, only: is_parallel
     19use mod_grid_phy_lmdz,  only: nbp_lon, nbp_lat
     20
     21implicit none
     22
     23! This routine computes slope inclination and orientation for the GCM (callslope=.true. in callphys.def)
     24! It works fine with a non-regular grid for zoomed simulations.
     25! slope inclination angle (deg)  0 == horizontal, 90 == vertical
     26! slope orientation angle (deg)  0 == Northward,  90 == Eastward, 180 == Southward, 270 == Westward
     27! TN 04/1013
     28
     29! Input arguments
     30integer,                intent(in) :: ngrid  ! nnumber of atmospheric columns
     31real, dimension(ngrid), intent(in) :: geopot ! geopotential on phy grid
     32
     33! Local variables
     34real, dimension(nbp_lon,nbp_lat) :: topogrid           ! topography on lat/lon grid with poles and only one -180/180 point
     35real, dimension(nbp_lon,nbp_lat) :: latigrid, longgrid ! meshgrid of latitude and longitude values (radians)
     36real, dimension(nbp_lon,nbp_lat) :: gradx              ! x: latitude-wise topography gradient,  increasing northward
     37real, dimension(nbp_lon,nbp_lat) :: grady              ! y: longitude-wise topography gradient, increasing westward
     38real                             :: theta_val          ! slope inclination
     39real                             :: psi_val            ! slope orientation
     40integer                          :: i, j, ig0
     41integer                          :: id2, idm1          ! a trick to compile testphys1d with debug option
     42
     43if (is_parallel) then
     44    ! This routine only works in serial mode so stop now.
     45    write(*,*) "getslopes Error: this routine is not designed to run in parallel"
     46    call abort_physic("getslopes",'cannot be run in parallel',1)
     47endif
     48
     49id2  = 2
     50idm1 = nbp_lon-1
     51
     52! rearrange topography on a 2d array
     53do j = 2,nbp_lat-1
     54    ig0 = 1 + (j - 2)*nbp_lon
     55    do i = 1,nbp_lon
     56        topogrid(i,j) = geopot(ig0 + i)/g
     57        latigrid(i,j) = latitude(ig0 + i)
     58        longgrid(i,j) = longitude(ig0 + i)
     59    enddo
     60enddo
     61
     62! poles:
     63topogrid(:,1) = geopot(1)/g
     64latigrid(:,1) = latitude(1)
     65longgrid(:,1) = longitude(1)
     66topogrid(:,nbp_lat) = geopot(ngrid)/g
     67latigrid(:,nbp_lat) = latitude(ngrid)
     68longgrid(:,nbp_lat) = longitude(ngrid)
     69
     70! compute topography gradient
     71! topogrid and rad are both in meters
     72do j = 2,nbp_lat - 1
     73    do i=1,nbp_lon
     74        gradx(i,j) = (topogrid(i,j + 1) - topogrid(i,j - 1))/(latigrid(i,j + 1)-latigrid(i,j - 1))
     75        gradx(i,j) = gradx(i,j)/rad
     76    enddo
     77    grady(1,j) = (topogrid(id2,j) - topogrid(nbp_lon,j))/(2*pi + longgrid(id2,j) - longgrid(nbp_lon,j))
     78    grady(1,j) = grady(1,j) / rad
     79    grady(nbp_lon,j) = (topogrid(1,j) - topogrid(idm1,j))/(2*pi + longgrid(1,j) - longgrid(idm1,j))
     80    grady(nbp_lon,j) = grady(nbp_lon,j)/rad
     81    do i = 2,nbp_lon - 1
     82        grady(i,j) = (topogrid(i + 1,j) - topogrid(i-1,j))/(longgrid(i + 1,j) - longgrid(i - 1,j))
     83        grady(i,j) = grady(i,j)/rad
     84    enddo
     85enddo
     86
     87! poles:
     88gradx(:,1) = 0.
     89grady(:,1) = 0.
     90gradx(:,nbp_lat) = 0.
     91grady(:,nbp_lat) = 0.
     92
     93! compute slope inclination and orientation:
     94theta_sl = 0.
     95psi_sl  = 0.
     96do j = 2,nbp_lat - 1
     97    do i = 1,nbp_lon
     98        ig0 = 1 + (j - 2)*nbp_lon
     99
     100        theta_val = atan(sqrt((gradx(i,j))**2 + (grady(i,j))**2))
     101
     102        psi_val = 0.
     103        if (gradx(i,j) /= 0.) psi_val = -pi/2. - atan(grady(i,j)/gradx(i,j))
     104        if (gradx(i,j) >= 0.) psi_val = psi_val - pi
     105        psi_val = 3*pi/2. - psi_val
     106        psi_val = psi_val*180./pi
     107        psi_val = modulo(psi_val,360.)
     108
     109        theta_sl(ig0 + i) = theta_val
     110        psi_sl(ig0 + i)   = psi_val
     111    enddo
     112enddo
     113
     114end subroutine getslopes
     115
     116!=======================================================================
     117
     118SUBROUTINE param_slope(csza,declin,rho,latitude,taudust,albedo,theta_s,psi_s,fdir_0,ftot_0,ftot)
     119!***********************************************************************
     120!
     121! SUBROUTINE:
     122! param_slope
     123!
     124!
     125! PURPOSE:
     126! computes total solar irradiance on a given Martian slope
     127!
     128!
     129! INPUTS:
     130! csza          cosine solar zenith angle
     131! declin        sun declination (rad)
     132! rho           sun right ascension (rad)
     133! latitude      latitude (deg)
     134! taudust       dust optical depth at reference wavelength 0.67 mic.
     135! albedo        spectrally integrated surface Lambertian reflection albedo
     136! theta_s       slope inclination angle (deg)
     137!                       0 is horizontal, 90 is vertical
     138! phi_s         slope azimuth (deg)
     139!                       0 >> Northward
     140!                       90 >> Eastward
     141!                       180 >> Southward
     142!                       270 >> Westward
     143! ftot_0        spectrally integrated total irradiance on an horizontal surface (W/m2)
     144! fdir_0        spectrally integrated direct irradiance on an horizontal surface (W/m2)
     145!
     146!
     147! OUTPUTS:
     148! ftot          spectrally integrated total irradiance on the slope (W/m2)
     149!
     150! REFERENCE:
     151! "Fast and accurate estimation of irradiance on Martian slopes"
     152! A. Spiga & F. Forget
     153! .....
     154!
     155! AUTHOR:
     156! A. Spiga (spiga@lmd.jussieu.fr)
     157! March 2008
     158!
     159!***********************************************************************
     160
     161use comcstfi_h, only: pi
     162
     163implicit none
     164
     165! Input arguments
     166real, intent(in) :: csza, declin, rho, latitude, taudust, theta_s, psi_s, albedo, ftot_0 , fdir_0
     167
     168! Output arguments
     169real, intent(out) :: ftot
     170
     171! Local variables
     172real                 :: deg2rad, a, mu_s, sigma_s, fdir, fscat, fscat_0, fref, ratio
     173real, dimension(4,2) :: mat_M, mat_N, mat_T
     174real, dimension(2)   :: g_vector
     175real, dimension(4)   :: s_vector
     176!***********************************************************************
     177! Prerequisite
     178deg2rad = pi/180.
     179if ((theta_s > 90.) .or. (theta_s < 0.)) then
     180    write(*,*) 'please set theta_s between 0 and 90', theta_s
     181    call abort_physic("param_slope","invalid theta_s",1)
     182endif
     183
     184! Solar Zenith angle (radian)
     185if (csza < 0.01) then
     186    !print *, 'sun below horizon'
     187    !fdir_0=0.
     188    fdir    = 0.
     189    fscat_0 = 0.
     190    fscat   = 0.   
     191    fref    = 0.
     192else
     193! Low incidence fix
     194!    if (csza < 0.15) csza = 0.15
     195
     196! 'Slope vs Sun' azimuth (radian)
     197    if (cos(declin)*sin(rho) == 0. .and. &
     198        sin(deg2rad*latitude)*cos(declin)*cos(rho) - cos(deg2rad*latitude)*sin(declin) == 0.) then
     199        a = deg2rad*psi_s ! some compilator need specfying value for atan2(0,0) 
     200    else
     201        a = deg2rad*psi_s + atan2(cos(declin)*sin(rho),sin(deg2rad*latitude)*cos(declin)*cos(rho)-cos(deg2rad*latitude)*sin(declin))
     202    endif
     203
     204! Cosine of slope-sun phase angle
     205    mu_s = csza*cos(deg2rad*theta_s) - cos(a)*sin(deg2rad*theta_s)*sqrt(1-csza**2)
     206    if (mu_s <= 0.) mu_s = 0.
     207
     208! Sky-view factor
     209    sigma_s=0.5*(1. + cos(deg2rad*theta_s))
     210
     211! Direct flux on the slope
     212    fdir = fdir_0*mu_s/csza
     213
     214! Reflected flux on the slope
     215    fref = albedo*(1 - sigma_s)*ftot_0
     216
     217! Scattered flux on a flat surface
     218    fscat_0 = ftot_0 - fdir_0
     219
     220! Scattering vector (slope vs sky)
     221    s_vector = (/ 1., exp(-taudust), sin(deg2rad*theta_s), sin(deg2rad*theta_s)*exp(-taudust) /)
     222
     223! Geometry vector (slope vs sun)
     224    g_vector = (/ mu_s/csza, 1. /)
     225
     226! Coupling matrix
     227    if (csza >= 0.5) then
     228        mat_M(:,1) = (/ -0.264,  1.309,  0.208, -0.828 /)
     229        mat_M(:,2) = (/  1.291*sigma_s, -1.371*sigma_s, -0.581,  1.641 /)
     230        mat_N(:,1) = (/  0.911, -0.777, -0.223,  0.623 /)
     231        mat_N(:,2) = (/ -0.933*sigma_s,  0.822*sigma_s,  0.514, -1.195 /)
     232    else
     233        mat_M(:,1) = (/ -0.373,  0.792, -0.095,  0.398 /)
     234        mat_M(:,2) = (/  1.389*sigma_s, -0.794*sigma_s, -0.325,  0.183 /)
     235        mat_N(:,1) = (/  1.079,  0.275,  0.419, -1.855 /)
     236        mat_N(:,2) = (/ -1.076*sigma_s, -0.357*sigma_s, -0.075,  1.844 /)
     237    endif
     238    mat_T = mat_M + csza*mat_N
     239
     240! Scattered flux slope ratio
     241    if (deg2rad*theta_s <= 0.0872664626) then
     242    ! low angles
     243        s_vector = (/ 1., exp(-taudust) , sin(0.0872664626), sin(0.0872664626)*exp(-taudust) /)
     244        ratio = dot_product(matmul(s_vector, mat_T),g_vector)
     245        ratio = 1. + (ratio - 1.)*deg2rad*theta_s/0.0872664626
     246    else
     247    ! general case
     248    ratio = dot_product(matmul(s_vector,mat_T),g_vector) 
     249    ! NB: ratio = dot_product(s_vector,matmul(mat_T,g_vector)) is equivalent
     250    endif
     251
     252! Scattered flux on the slope
     253    fscat = ratio*fscat_0
     254endif ! if (csza < 0.01)
     255
     256! Total flux on the slope
     257ftot = fdir + fref + fscat
     258
     259! Display results
     260!  print *, 'sca component 0 ', ftot_0-fdir_0
     261!  print *, 'dir component 0 ', fdir_0
     262!  print *, 'scattered component ', fscat
     263!  print *, 'direct component ', fdir
     264!  print *, 'reflected component ', fref
     265
     266END SUBROUTINE param_slope
     267
     268!=======================================================================
     269
     270SUBROUTINE ini_slope_mod(ngrid)
     271
     272implicit none
     273
     274integer, intent(in) :: ngrid ! number of atmospheric columns
     275
     276allocate(theta_sl(ngrid))
     277allocate(psi_sl(ngrid))
     278
     279END SUBROUTINE ini_slope_mod
     280
     281!=======================================================================
     282
     283SUBROUTINE end_slope_mod
     284
     285implicit none
     286
     287if (allocated(theta_sl)) deallocate(theta_sl)
     288if (allocated(psi_sl)) deallocate(psi_sl)
     289
     290END SUBROUTINE end_slope_mod
     291
     292END MODULE slope_mod
Note: See TracChangeset for help on using the changeset viewer.