Changeset 4035 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Jan 30, 2026, 12:40:55 PM (4 weeks ago)
Author:
emillour
Message:

Mars PCM:
Cleanup and correction around thermospheric processes: ensure that the
tendencies from previous processes are included before being sent to
the next process (was not the case for conduction).
EM

Location:
trunk/LMDZ.MARS/libf
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90

    r4023 r4035  
    66
    77      SUBROUTINE euvheat(ngrid,nlayer,nq,pt,pdt,pplev,pplay,zzlay, &
    8            mu0,ptimestep,ptime,zday,pq,pdq,pdteuv)
     8           mu0,ptimestep,zday,pq,pdq,pdteuv)
    99
    1010      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,         &
     
    4747      integer,intent(in) :: nlayer ! number of atmospheric layers
    4848      integer,intent(in) :: nq ! number of advected tracers
    49       real :: pt(ngrid,nlayer)
    50       real :: pdt(ngrid,nlayer)
    51       real :: pplev(ngrid,nlayer+1)
    52       real :: pplay(ngrid,nlayer)
    53       real :: zzlay(ngrid,nlayer)
    54       real :: mu0(ngrid)
    55       real :: ptimestep,ptime
    56       real :: zday
    57       real :: pq(ngrid,nlayer,nq)
    58       real :: pdq(ngrid,nlayer,nq)
    59 
    60       real :: pdteuv(ngrid,nlayer)
     49      real,intent(in) :: pt(ngrid,nlayer) ! input temperature (K)
     50      real,intent(in) :: pdt(ngrid,nlayer) ! input tendency on temperatur (K/s)
     51      real,intent(in) :: pplev(ngrid,nlayer+1) ! pressure (Pa) at layer interfaces
     52      real,intent(in) :: pplay(ngrid,nlayer) ! pressure (Pa) at mid-layer
     53      real,intent(in) :: zzlay(ngrid,nlayer) ! mid-layer altitude (m)
     54      real,intent(in) :: mu0(ngrid) ! cosine of solar zenith angle
     55      real,intent(in) :: ptimestep ! physics time step (s)
     56      real,intent(in) :: zday ! date (sols since Ls=0)
     57      real,intent(in) :: pq(ngrid,nlayer,nq) ! input tracer mmr (kg/kg_air)
     58      real,intent(in) :: pdq(ngrid,nlayer,nq) ! input tendencies on tracers (kg/kg_air/s)
     59
     60      real,intent(out) :: pdteuv(ngrid,nlayer) ! temperature tendency (K/s)
    6161!
    6262!    Local variables :
  • trunk/LMDZ.MARS/libf/aeronomars/moldiff_MPF.F90

    r4021 r4035  
    88
    99subroutine moldiff_MPF(ngrid,nlayer,nq,pplay,pplev,pt,pdt,pq,pdq,&
    10                        ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff,&
     10                       ptimestep,pdqdiff,&
    1111                       PhiEscH,PhiEscH2,PhiEscD)
    1212
     
    3232      integer,intent(in) :: nlayer ! number of atmospheric layers
    3333      integer,intent(in) :: nq ! number of advected tracers
    34       real ptimestep
    35       real pplay(ngrid,nlayer)
    36       real zzlay(ngrid,nlayer)
    37       real pplev(ngrid,nlayer+1)
    38       real pq(ngrid,nlayer,nq)
    39       real pdq(ngrid,nlayer,nq)
    40       real pt(ngrid,nlayer)
    41       real pdt(ngrid,nlayer)
    42       real pdteuv(ngrid,nlayer)
    43       real pdtconduc(ngrid,nlayer)
    44       real pdqdiff(ngrid,nlayer,nq)
    45       real*8 PhiEscH,PhiEscH2,PhiEscD
     34      real,intent(in) :: ptimestep ! physics time step (s)
     35      real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
     36      real,intent(in) :: pplev(ngrid,nlayer+1) ! pressure (Pa) at layer boundaries
     37      real,intent(in) :: pq(ngrid,nlayer,nq) ! input tracer mmr (kg/kg_air)
     38      real,intent(in) :: pdq(ngrid,nlayer,nq) ! input tendency on tracers (kg/kg_air/s)
     39      real,intent(in) :: pt(ngrid,nlayer) ! input temperature (K)
     40      real,intent(in) :: pdt(ngrid,nlayer) ! input tendency on temperature (K/s)
     41      real,intent(out) :: pdqdiff(ngrid,nlayer,nq) ! tendency on tracers (kg/kg_air/s)
     42      real*8,intent(out) :: PhiEscH ! escape flux of H (.../s)
     43      real*8,intent(out) :: PhiEscH2 ! escape flux of H2 (.../s)
     44      real*8,intent(out) :: PhiEscD ! escape flux of D (.../s)
    4645!
    4746! Local
     
    271270!     &  ,tt,ptimestep,nlayer,ig)
    272271        do l=1,nlayer
    273           tt(l)=pt(ig,l)*1D0+(pdt(ig,l)*dble(ptimestep)+                &
    274                               pdtconduc(ig,l)*dble(ptimestep)+          &
    275                               pdteuv(ig,l)*dble(ptimestep))
     272          tt(l)=pt(ig,l)*1D0+(pdt(ig,l)*dble(ptimestep))
    276273          ! to cach Nans...
    277274          if (tt(l).ne.tt(l)) then
    278275            print*,'Err TMNEW',ig,l,tt(l),pt(ig,l),                     &
    279               pdt(ig,l),pdtconduc(ig,l),pdteuv(ig,l),dble(ptimestep)
     276              pdt(ig,l),dble(ptimestep)
    280277          endif
    281278        enddo ! of do l=1,nlayer
     
    338335        Print*,'Zmax too high',ig,zmax,zmin
    339336        do l=1,nlayer
    340         print*,'old',zz(l),pt(ig,l),pdteuv(ig,l),pdq(ig,l,:)
     337        print*,'old',zz(l),pt(ig,l),pdt(ig,l),pdq(ig,l,:)
    341338        print*,'l',l,rhot(l),tt(l),pp(l),massemoy(l),qq(l,:)
    342339        enddo
  • trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90

    r4018 r4035  
    88
    99subroutine moldiff_red(ngrid,nlayer,nq,pplay,pplev,pt,pdt,pq,pdq,&
    10                        ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff,&
     10                       ptimestep,pdqdiff,&
    1111                       PhiEscH,PhiEscH2,PhiEscD)
    1212
     
    3131      integer,intent(in) :: nlayer ! number of atmospheric layers
    3232      integer,intent(in) :: nq ! number of advected tracers
    33       real ptimestep
    34       real pplay(ngrid,nlayer)
    35       real zzlay(ngrid,nlayer)
    36       real pplev(ngrid,nlayer+1)
    37       real pq(ngrid,nlayer,nq)
    38       real pdq(ngrid,nlayer,nq)
    39       real pt(ngrid,nlayer)
    40       real pdt(ngrid,nlayer)
    41       real pdteuv(ngrid,nlayer)
    42       real pdtconduc(ngrid,nlayer)
    43       real pdqdiff(ngrid,nlayer,nq)
    44       real*8 PhiEscH,PhiEscH2,PhiEscD
     33      real,intent(in) :: ptimestep ! physics time step (s)
     34      real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
     35      real,intent(in) :: pplev(ngrid,nlayer+1) ! pressure (Pa) at layer boundaries
     36      real,intent(in) :: pq(ngrid,nlayer,nq) ! input tracer mmr (kg/kg_air)
     37      real,intent(in) :: pdq(ngrid,nlayer,nq) ! input tendency on tracers (kg/kg_air/s)
     38      real,intent(in) :: pt(ngrid,nlayer) ! input temperature (K)
     39      real,intent(in) :: pdt(ngrid,nlayer) ! input tendency on temperature (K/s)
     40      real,intent(out) :: pdqdiff(ngrid,nlayer,nq) ! tendency on tracers (kg/kg_air/s)
     41      real*8,intent(out) :: PhiEscH,PhiEscH2,PhiEscD
    4542!
    4643! Local
     
    255252!     &  ,tt,ptimestep,nlayer,ig)
    256253        do l=1,nlayer
    257           tt(l)=pt(ig,l)*1D0+(pdt(ig,l)*dble(ptimestep)+                &
    258                               pdtconduc(ig,l)*dble(ptimestep)+          &
    259                               pdteuv(ig,l)*dble(ptimestep))
     254          tt(l)=pt(ig,l)*1D0+(pdt(ig,l)*dble(ptimestep))
    260255          ! to cach Nans...
    261256          if (tt(l).ne.tt(l)) then
    262257            print*,'Err TMNEW',ig,l,tt(l),pt(ig,l),                     &
    263               pdt(ig,l),pdtconduc(ig,l),pdteuv(ig,l),dble(ptimestep)
     258              pdt(ig,l),dble(ptimestep)
    264259          endif
    265260        enddo ! of do l=1,nlayer
     
    322317        Print*,'Zmax too high',ig,zmax,zmin
    323318        do l=1,nlayer
    324         print*,'old',zz(l),pt(ig,l),pdteuv(ig,l),pdq(ig,l,:)
     319        print*,'old',zz(l),pt(ig,l),pdt(ig,l),pdq(ig,l,:)
    325320        print*,'l',l,rhot(l),tt(l),pp(l),massemoy(l),qq(l,:)
    326321        enddo
  • trunk/LMDZ.MARS/libf/aeronomars/molvis.F

    r3466 r4035  
    66
    77      SUBROUTINE molvis(ngrid,nlayer,ptimestep,
    8      &            pplay,pplev,pt,pdteuv,pdtconduc
    9      $           ,pvel,tsurf,zzlev,zzlay,zdvelmolvis)
     8     &            pplay,pplev,pt,pdt,
     9     $            pvel,tsurf,zzlev,zzlay,zdvelmolvis)
    1010     
    1111      use conc_mod, only: cpnew, Akknew, rnew
     
    3131      integer,intent(in) :: ngrid ! number of atmospheric columns
    3232      integer,intent(in) :: nlayer ! number of atmospheric layers
    33       REAL ptimestep
    34       REAL pplay(ngrid,nlayer)
    35       REAL pplev(ngrid,nlayer+1)
    36       REAL zzlay(ngrid,nlayer)
    37       REAL zzlev(ngrid,nlayer+1)
    38       real pt(ngrid,nlayer)
    39       real tsurf(ngrid)
    40       REAL pvel(ngrid,nlayer)
    41       REAL pdvel(ngrid,nlayer)
    42       real pdteuv(ngrid,nlayer)
    43       real pdtconduc(ngrid,nlayer)
     33      REAL,INTENT(IN) :: ptimestep ! physics time step (s)
     34      REAL,INTENT(IN) :: pplay(ngrid,nlayer)
     35      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
     36      REAL,INTENT(IN) :: zzlay(ngrid,nlayer)
     37      REAL,INTENT(IN) :: zzlev(ngrid,nlayer+1)
     38      real,intent(in) :: pt(ngrid,nlayer) ! input temperature (K)
     39      real,intent(in) :: pdt(ngrid,nlayer) ! input tendency on temperature (K/s)
     40      real,intent(in) :: tsurf(ngrid) ! input ssurface temperature (K)
     41      REAL,INTENT(IN) :: pvel(ngrid,nlayer) ! wind velocity (m/s)
    4442
    45       real zdvelmolvis(ngrid,nlayer)
     43      real,intent(out) :: zdvelmolvis(ngrid,nlayer) ! tendency on velocity (m/s/s)
    4644
    4745c   local:
     
    9997      do ig=1,ngrid
    10098
    101         zt(1)=pt(ig,1)+(pdteuv(ig,1)+pdtconduc(ig,1))*ptimestep
     99        zt(1)=pt(ig,1)+(pdt(ig,1))*ptimestep
    102100        zvel(1)=pvel(ig,1)
    103101        zlay(1)=zzlay(ig,1)
     
    105103
    106104        do l=2,nz
    107           zt(l)=pt(ig,l)+(pdteuv(ig,l)+pdtconduc(ig,l))*ptimestep
     105          zt(l)=pt(ig,l)+(pdt(ig,l))*ptimestep
    108106          zvel(l)=pvel(ig,l)
    109107          zlay(l)=zzlay(ig,l)
  • trunk/LMDZ.MARS/libf/aeronomars/thermosphere_mod.F

    r4024 r4035  
    1010      subroutine thermosphere(ngrid,nlayer,nq,
    1111     &     pplev,pplay,dist_sol,
    12      $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
    13      &     pt,pq,pu,pv,pdt,pdq,
     12     $     mu0,ptimestep,zday,tsurf,zzlev,zzlay,
     13     &     pt,pq,pu,pv,pdt,pdq,pdu,pdv,
    1414     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
    1515     $     PhiEscH,PhiEscH2,PhiEscD)
     
    3131      integer,intent(in) :: nlayer ! number of atmospheric layers
    3232      integer,intent(in) :: nq ! number of advected tracers
    33       REAL,INTENT(in) :: pplay(ngrid,nlayer)
    34       REAL,INTENT(in) :: pplev(ngrid,nlayer+1)
    35       REAL,INTENT(in) :: zzlay(ngrid,nlayer)
    36       REAL,INTENT(in) :: zzlev(ngrid,nlayer+1)
    37       REAL,INTENT(in) :: pt(ngrid,nlayer)
    38       REAL,INTENT(in) :: zday
    39       REAL,INTENT(in) :: dist_sol
    40       REAL,INTENT(in) :: mu0(ngrid)
    41       REAL,INTENT(in) :: pq(ngrid,nlayer,nq)
    42       REAL,INTENT(in) :: ptimestep
    43       REAL,INTENT(in) :: ptime
    44       REAL,INTENT(in) :: tsurf(ngrid)
    45       REAL,INTENT(in) :: pu(ngrid,nlayer),pv(ngrid,nlayer)
    46       REAL,INTENT(in) :: pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
     33      REAL,INTENT(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
     34      REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
     35      REAL,INTENT(in) :: zzlay(ngrid,nlayer) ! altitude of mid-layer (m)
     36      REAL,INTENT(in) :: zzlev(ngrid,nlayer+1) ! altitude of layer boundaries (m)
     37      REAL,INTENT(in) :: pt(ngrid,nlayer) ! input temperature (K)
     38      REAL,INTENT(in) :: zday ! Mars date (sols since Ls=0)
     39      REAL,INTENT(in) :: dist_sol ! Sun-Mars distance (AU)
     40      REAL,INTENT(in) :: mu0(ngrid) ! cosine of solar zenith angle
     41      REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! input tracer mixing ratio (kg/kg_air)
     42      REAL,INTENT(in) :: ptimestep ! physics time step (s)
     43      REAL,INTENT(in) :: tsurf(ngrid) ! surface temperature (K)
     44      REAL,INTENT(in) :: pu(ngrid,nlayer) ! input zonal wind (m/s)
     45      REAL,INTENT(in) :: pv(ngrid,nlayer) ! input meridional wind (m/s)
    4746
    48       REAL,INTENT(out) :: zdteuv(ngrid,nlayer)
    49       REAL,INTENT(out) :: zdtconduc(ngrid,nlayer)
    50       REAL,INTENT(out) :: zdumolvis(ngrid,nlayer)
    51       REAL,INTENT(out) :: zdvmolvis(ngrid,nlayer)
    52       REAL,INTENT(out) :: zdqmoldiff(ngrid,nlayer,nq)
    53       REAL*8,INTENT(out) :: PhiEscH,PhiEscH2,PhiEscD
     47      REAL,INTENT(out) :: zdteuv(ngrid,nlayer) ! tendency due to EUV (K/s)
     48      REAL,INTENT(out) :: zdtconduc(ngrid,nlayer) ! tendency due to conduction (K/s)
     49      REAL,INTENT(out) :: zdumolvis(ngrid,nlayer) ! tendency due to viscosity (m/s/s)
     50      REAL,INTENT(out) :: zdvmolvis(ngrid,nlayer) ! tendency due to viscosity (m/s/s)
     51      REAL,INTENT(out) :: zdqmoldiff(ngrid,nlayer,nq) ! tendency due to diffusion (kg/kg_air/s)
     52      REAL*8,INTENT(out) :: PhiEscH ! total escape flux of H (s-1)
     53      REAL*8,INTENT(out) :: PhiEscH2 ! total escape flux of H2 (s-1)
     54      REAL*8,INTENT(out) :: PhiEscD ! total escape flux of D (s-1)
     55
     56      REAL,INTENT(inout) :: pdt(ngrid,nlayer) ! tendency on temperature (K/s)
     57      REAL,INTENT(inout) :: pdq(ngrid,nlayer,nq) ! tendency on tracers (kg/kg_air/s)
     58      REAL,INTENT(inout) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s/s)
     59      REAL,INTENT(inout) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s/s)
    5460
    5561      INTEGER :: l,ig
     
    6874
    6975      ! initialize tendencies to zero in all cases
    70       ! (tendencies are added later on, even if parametrization is not called)
    7176      zdteuv(1:ngrid,1:nlayer)=0
    7277      zdtconduc(1:ngrid,1:nlayer)=0
     
    7782      if (calleuv) then
    7883        call euvheat(ngrid,nlayer,nq,pt,pdt,pplev,pplay,zzlay,
    79      $               mu0,ptimestep,ptime,zday,pq,pdq,zdteuv)
     84     $               mu0,ptimestep,zday,pq,pdq,zdteuv)
     85        ! update tendency on temperature
     86        pdt(:,:)=pdt(:,:)+zdteuv(:,:)
    8087      endif
    8188
    8289      if (callconduct) THEN
    83         call conduction(ngrid,nlayer,ptimestep,pplay,pplev,pt,zdteuv,
     90        call conduction(ngrid,nlayer,ptimestep,pplay,pplev,pt,pdt,
    8491     $                   tsurf,zzlev,zzlay,zdtconduc)
     92        ! update tendency on temperature
     93        pdt(:,:)=pdt(:,:)+zdtconduc(:,:)
    8594      endif
    8695
    8796      if (callmolvis) THEN
    8897        call molvis(ngrid,nlayer,ptimestep,pplay,pplev,pt,
    89      &                zdteuv,zdtconduc,pu,
    90      $                   tsurf,zzlev,zzlay,zdumolvis)
     98     &                pdt,pu,tsurf,zzlev,zzlay,zdumolvis)
    9199        call molvis(ngrid,nlayer,ptimestep,pplay,pplev,pt,
    92      &                zdteuv,zdtconduc,pv,
    93      $                   tsurf,zzlev,zzlay,zdvmolvis)
     100     &                pdt,pv,tsurf,zzlev,zzlay,zdvmolvis)
     101        ! update tendencies on winds
     102        pdu(:,:)=pdu(:,:)+zdumolvis(:,:)
     103        pdv(:,:)=pdv(:,:)+zdvmolvis(:,:)
    94104      endif
    95105
     
    99109        call moldiff_red(ngrid,nlayer,nq,
    100110     &                   pplay,pplev,pt,pdt,pq,pdq,ptimestep,
    101      &                   zzlay,zdteuv,zdtconduc,zdqmoldiff,
    102      &                   PhiEscH,PhiEscH2,PhiEscD)
     111     &                   zdqmoldiff,PhiEscH,PhiEscH2,PhiEscD)
    103112       else
    104113        ! new MPF (modified pass flow) scheme
    105114        call moldiff_MPF(ngrid,nlayer,nq,
    106115     &                   pplay,pplev,pt,pdt,pq,pdq,ptimestep,
    107      &                   zzlay,zdteuv,zdtconduc,zdqmoldiff,
    108      &                   PhiEscH,PhiEscH2,PhiEscD)
     116     &                   zdqmoldiff,PhiEscH,PhiEscH2,PhiEscD)
    109117       endif ! of if (moldiff_scheme==1)
    110       endif
     118
     119       ! update tendencies on tracers
     120       pdq(:,:,:)=pdq(:,:,:)+zdqmoldiff(:,:,:)
     121       
     122      endif ! of if (callmoldiff)
    111123
    112124      end subroutine thermosphere
  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r4029 r4035  
    21302130c     ------------------
    21312131
    2132 !#ifndef MESOSCALE
    21332132c        --------------
    21342133c        photochemistry :
     
    22002199
    22012200         END IF  ! of IF (photochem)
    2202 !#endif
    2203 
    2204 
    2205 !#ifndef MESOSCALE
     2201
     2202
    22062203c-----------------------------------------------------------------------
    22072204c   10. THERMOSPHERE CALCULATION
     
    22102207      if (callthermos) then
    22112208        call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
    2212      $     mu0,ptimestep,ptime,zday,tsurf_meshavg,zzlev,zzlay,
    2213      &     pt,pq,pu,pv,pdt,pdq,
    2214      $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
    2215      $     PhiEscH,PhiEscH2,PhiEscD)
    2216 
    2217         DO l=1,nlayer
    2218           DO ig=1,ngrid
    2219             dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
    2220             pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l)
    2221             pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
    2222             pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
    2223             DO iq=1, nq
    2224               pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
    2225             ENDDO
    2226           ENDDO
    2227         ENDDO
     2209     &     mu0,ptimestep,zday,tsurf_meshavg,zzlev,zzlay,
     2210     &     pt,pq,pu,pv,pdt,pdq,pdu,pdv,
     2211     &     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
     2212     &     PhiEscH,PhiEscH2,PhiEscD)
     2213       
     2214        ! tendencies zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff
     2215        ! already included in pdt,pdu,pdv,pdq by thermosphere
     2216       
     2217        dtrad(1:ngrid,1:nlayer)=dtrad(1:ngrid,1:nlayer)
     2218     &                               +zdteuv(1:ngrid,1:nlayer)
    22282219
    22292220      endif ! of if (callthermos)
    2230 !#endif
    22312221
    22322222c-----------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.