Changeset 3115 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Nov 3, 2023, 6:02:56 PM (13 months ago)
Author:
llange
Message:

Mars PCM

  • Introducing the possibily to compute water adsorption / desorption -routine soilwater.F90) -. FOR NOW IT WORKS ONLY IN 1D, DON'T TEST IN 3D (by default, adsorption is not called, if not specified in the callphys.def). By default, when using the adsorption, the number of subtimestep used in the water subimation scheme is 1 (otherwise, it crashes)
  • New grid for the soil layers (better resolution) as it is needed to solve the diffusion equations. It does not increase the CPU time.

PYM & LL

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

Legend:

Unmodified
Added
Removed

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

    r3113 r3115  
    33implicit none
    44! nsoilmx : number of subterranean layers
    5   integer, parameter :: nsoilmx = 18
     5  integer, parameter :: nsoilmx = 57
    66
    77  real,save,allocatable,dimension(:) :: layer      ! soil layer depths

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

    r3113 r3115  
    670670! -----------------
    671671do isoil = 0,nsoil - 1
    672     mlayer(isoil) = 2.e-4*(2.**(isoil - 0.5)) ! mid-layer depth
    673     layer(isoil + 1) = 2.e-4*(2.**isoil)      ! layer depth
     672    mlayer(isoil) = 2.e-4*(1+isoil**2.9*(1-exp(-real(isoil)/20.))) ! mid layer depth
    674673enddo
     674
     675do isoil = 1,nsoil - 1
     676    layer(isoil)=(mlayer(isoil)+mlayer(isoil-1))/2
     677enddo
     678layer(nsoil)=2*mlayer(nsoil-1) - mlayer(nsoil-2)
    675679
    676680! Initialize traceurs

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

    r3113 r3115  
    15071507         zcdv(:) = 0.
    15081508
    1509          CALL vdifc(ngrid,nlayer,nq,zpopsk,
     1509         CALL vdifc(ngrid,nlayer,nsoilmx,nq,nqsoil,zpopsk,
    15101510     $        ptimestep,capcal,lwrite,
    15111511     $        zplay,zplev,zzlay,zzlev,z0,
    1512      $        pu,pv,zh,pq,tsurf,emis,qsurf,
     1512     $        pu,pv,zh,pq,tsurf,tsoil,emis,qsurf,qsoil,
    15131513     $        zdum1,zdum2,zdh,pdq,zflubid,
    15141514     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
    15151515     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
    15161516     &        zcondicea_co2microp,sensibFlux,
    1517      &        dustliftday,local_time,watercap,dwatercap_dif,
    1518      &        tsoil,nsoilmx)
     1517     &        dustliftday,local_time,watercap,dwatercap_dif)
    15191518
    15201519          DO ig=1,ngrid

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

    r3113 r3115  
    148148! 1.4 Build mlayer(), if necessary
    149149      if (interpol) then
    150       ! default mlayer distribution, following a power law:
    151       !  mlayer(k)=lay1*alpha**(k-1/2)
    152          lay1=2.e-4
    153          alpha=2
    154          do iloop=0,nsoil-1
    155            mlayer(iloop)=lay1*(alpha**(iloop-0.5))
    156          enddo
     150      !  mlayer(k)=lay1*(1+k^(2.9)*(1-exp(-k/20)), PYM grid
     151        lay1=2.e-4
     152        do iloop=0,nsoil-1
     153          mlayer(iloop) = lay1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
     154        enddo
    157155      endif
    158156
    159 ! 1.5 Build layer(); following the same law as mlayer()
    160       ! Assuming layer distribution follows mid-layer law:
    161       ! layer(k)=lay1*alpha**(k-1)
    162       lay1=sqrt(mlayer(0)*mlayer(1))
    163       alpha=mlayer(1)/mlayer(0)
    164       do iloop=1,nsoil
    165         layer(iloop)=lay1*(alpha**(iloop-1))
     157! 1.5 Build layer();
     158      do iloop=1,nsoil-1
     159        layer(iloop)=(mlayer(iloop)+mlayer(iloop-1))/2
    166160      enddo
    167 
     161      layer(nsoil)=2*mlayer(nsoil-1) - mlayer(nsoil-2)
    168162
    169163! 2. Volumetric heat capacity (note: it is declared in comsoil_h)

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

    r3111 r3115  
    55      CONTAINS
    66     
    7       SUBROUTINE vdifc(ngrid,nlay,nq,ppopsk,
     7      SUBROUTINE vdifc(ngrid,nlay,nsoil,nq,nqsoil,ppopsk,
    88     $                ptimestep,pcapcal,lecrit,
    99     $                pplay,pplev,pzlay,pzlev,pz0,
    10      $                pu,pv,ph,pq,ptsrf,pemis,pqsurf,
     10     $                pu,pv,ph,pq,ptsrf,ptsoil,pemis,pqsurf,qsoil,
    1111     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
    1212     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
    1313     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
    1414     $                hfmax,pcondicea_co2microp,sensibFlux,
    15      $                dustliftday,local_time,watercap, dwatercap_dif,
    16      $                tsoil,nsoil)
     15     $                dustliftday,local_time,watercap, dwatercap_dif)
     16
    1717      use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
    1818     &                      igcm_dust_submicron, igcm_h2o_vap,
    … …  
    3434      use microphys_h, only: To
    3535      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
    36       use comsoil_h, only: layer, mlayer
     36      use comsoil_h, only: layer, mlayer,adsorption_soil
    3737
    3838      IMPLICIT NONE
    … …  
    6464c   ----------
    6565
    66       INTEGER,INTENT(IN) :: ngrid,nlay,nsoil
     66      INTEGER,INTENT(IN) :: ngrid,nlay,nsoil,nqsoil
    6767      REAL,INTENT(IN) :: ptimestep
    6868      REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    … …  
    7878      REAL,INTENT(IN) :: pcapcal(ngrid,nslope)
    7979      REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
    80       REAL,INTENT(IN) :: tsoil(ngrid,nsoil,nslope)
     80      REAL,INTENT(IN) :: ptsoil(ngrid,nsoil,nslope)
    8181
    8282c    Argument added for condensation:
    … …  
    9999      real,intent(out) :: pdqsdif(ngrid,nq,nslope)
    100100      REAL,INTENT(in) :: dustliftday(ngrid)
     101      REAL,INTENT(inout) :: qsoil(ngrid,nsoil,nqsoil,nslope) !subsurface tracers
    101102      REAL,INTENT(in) :: local_time(ngrid)
    102103     
    … …  
    147148c     Subtimestep & implicit treatment of water vapor
    148149c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    149       REAL zdqsdif(ngrid) ! subtimestep pdqsdif for water ice
     150      REAL zdqsdif_surf(ngrid) ! subtimestep pdqsdif for water ice
    150151      REAL ztsrf(ngrid) ! temporary surface temperature in tsub
    151152      REAL zdtsrf(ngrid,nslope) ! surface temperature tendancy in tsub
    … …  
    227228
    228229      character*2 str2
     230
     231!! Subsurface exchanges
     232      LOGICAL :: exchange             ! boolean to check if exchange between the subsurface and the atmosphere can occurs
     233      REAL :: zdqsdif_regolith(ngrid,nslope) ! Flux from subsurface (positive pointing outwards) (kg/m^2/s)
     234      REAL zq1temp_regolith(ngrid)    ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg)
     235      REAL zdqsdif_tot(ngrid)         ! subtimestep pdqsdif for water ice
     236      LOGICAL :: writeoutput          ! boolean to say to soilexchange.F if we are at the last iteration and thus if he  can write in the diagsoil
    229237
    230238c    ** un petit test de coherence
    … …  
    283291      pdqsdif(1:ngrid,1:nq,1:nslope)=0
    284292      pdqsdif_tmp(1:ngrid,1:nq)=0
    285       zdqsdif(1:ngrid)=0
     293      zdqsdif_surf(1:ngrid)=0
    286294      dwatercap_dif(1:ngrid,1:nslope)=0
     295      zdqsdif_regolith(1:ngrid,1:nslope)=0
     296      zq1temp_regolith(1:ngrid)=0
     297      zdqsdif_tot(1:ngrid)=0
    287298!      h2o_ice_depth(1:ngrid,1:nslope)=5
    288299c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
    … …  
    944955c          make_tsub : sous pas de temps adaptatif
    945956c          la subroutine est a la fin du fichier
    946 
    947            call make_tsub(ngrid,pdtsrf(:,islope),zqsurf,
     957           if (adsorption_soil) then 
     958              nsubtimestep(:) = 1
     959           else
     960              call make_tsub(ngrid,pdtsrf(:,islope),zqsurf,
    948961     &                    ptimestep,dtmax,watercaptag,
    949962     &                    nsubtimestep)
    950 
     963           endif
    951964c           Calculation for turbulent exchange with the surface (for ice)
    952965c           initialization of ztsrf, which is surface temperature in
    … …  
    959972c           Debut du sous pas de temps
    960973            DO tsub=1,nsubtimestep(ig)
    961 
     974              if(tsub.eq.nsubtimestep(ig)) writeoutput = .true.
    962975c           C'est parti !
    963976             zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
    … …  
    9921005             zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
    9931006             if(old_wsublimation_scheme) then
    994                 zdqsdif(ig)=rho(ig)*dryness(ig)*zcdv(ig)
     1007                zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdv(ig)
    9951008     &                           *(zq1temp(ig)-qsat(ig))
    9961009             else
    997                 zdqsdif(ig)=rho(ig)*dryness(ig)*zcdh(ig)
     1010                zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdh(ig)
    9981011     &                           *(zq1temp(ig)-qsat(ig))
    9991012             endif
    1000 c             write(*,*)'subliming more than available frost:  qsurf!'
    1001 
     1013
     1014
     1015!!! Subsurface exchange
     1016! Check for subsurface exchanges
    10021017             if(.not.watercaptag(ig)) then
    1003                if ((-zdqsdif(ig)*subtimestep)
     1018                if (((-(zdqsdif_surf(ig))*
     1019     &             subtimestep).gt.zqsurf(ig))
     1020     &             .and.(pqsurf(ig,igcm_co2,islope).eq.0.)) then
     1021                      exchange = .true.
     1022                 else
     1023                      exchange = .false.
     1024                 endif
     1025             else
     1026                  exchange = .false.
     1027             endif
     1028                  zdqsdif_tot(ig) = zdqsdif_surf(ig)
     1029
     1030
     1031             if (adsorption_soil) then 
     1032
     1033                 call soilwater(1,nlay,nq,nsoil, nqsoil,
     1034     &                     ztsrf(ig),ptsoil(ig,:,islope),subtimestep,
     1035     &                     exchange,qsat(ig),zq_tmp_vap(ig,:,:),
     1036     &                     za(ig,:),zb(ig,:),zc(ig,:),zd(ig,:),
     1037     &                     zdqsdif_surf(ig), zqsurf(ig),
     1038     &                     qsoil(ig,:,:,islope), pplev(ig,1), rho(ig),
     1039     &                     writeoutput,zdqsdif_regolith(ig,islope),
     1040     &                     zq1temp_regolith(ig))
     1041
     1042
     1043
     1044                 if(.not.watercaptag(ig)) then
     1045                     if (exchange) then
     1046                         zq1temp(ig) = zq1temp_regolith(ig)
     1047                         zdqsdif_tot(ig)=
     1048     &                        -zqsurf(ig)/subtimestep
     1049                        zdqsdif_tot(ig) = zdqsdif_surf(ig)
     1050                     else
     1051                         zdqsdif_tot(ig) = zdqsdif_surf(ig) +
     1052     &                           zdqsdif_regolith(ig,islope) ! boundary condition = qsat, but pdqsdif is calculated to update qsurf (including loss of surface ice to the subsurface)
     1053                     endif  ! of "if exchange = true"
     1054                 endif  ! of "if not.watercaptag"
     1055              endif ! adsorption
     1056
     1057             if(.not.watercaptag(ig)) then
     1058               if ((-zdqsdif_tot(ig)*subtimestep)
    10041059     &            .gt.(zqsurf(ig))) then
    10051060               
    … …  
    10081063                IF(h2o_ice_depth(ig,islope) .gt. 0 .and. lag_layer)
    10091064     &           then
    1010                 zdqsdif(ig)=
     1065                zdqsdif_tot(ig)=
    10111066     &                        -zqsurf(ig)/subtimestep
    10121067                zqsurf(ig)=0
    … …  
    10271082     &                               +dist_down(ig,islope)
    10281083                 Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
    1029      &                          *tsoil(ig,lice-1,islope)
     1084     &                          *ptsoil(ig,lice-1,islope)
    10301085     &           /dist_sum(ig,islope))+
    1031      &           (dist_down(ig,islope)*tsoil(ig,lice,islope)
     1086     &           (dist_down(ig,islope)*ptsoil(ig,lice,islope)
    10321087     &           /dist_sum(ig,islope))
    10331088                        ELSE
    1034                                 Tice(ig,islope)=tsoil(ig,1,islope)
     1089                                Tice(ig,islope)=ptsoil(ig,1,islope)
    10351090                        ENDIF
    10361091              call watersat(1,Tice(ig,1),pplev(ig,1)
    … …  
    10791134             !write(*,*) "zdq_ssi*t", zdqsdif_ssi(ig)*subtimestep
    10801135             if (zdqsdif_ssi(ig,1)<0) then
    1081                !zdqsdif(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep)
    1082                zdqsdif(ig)=zdqsdif(ig)+zdqsdif_ssi(ig,1)
     1136               !zdqsdif_surf(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep)
     1137               zdqsdif_tot(ig)=zdqsdif_tot(ig)+zdqsdif_ssi(ig,1)
    10831138                call write_output('zdq_zdqssi',
    1084      &                '','',zdqsdif(ig)+zdqsdif_ssi(ig,1))
     1139     &                '','',zdqsdif_tot(ig)+zdqsdif_ssi(ig,1))
    10851140             endif
    10861141              ELSE
    … …  
    10881143c             pdqsdif < 0 : ice subliming
    10891144c             write(*,*) "subliming more than available frost:  qsurf!"
    1090                   zdqsdif(ig)=
     1145                  zdqsdif_tot(ig)=
    10911146     &                        -zqsurf(ig)/subtimestep
    10921147c                write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
    … …  
    10941149                 zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,igcm_h2o_vap)+
    10951150     $           zb(ig,2)*zc(ig,2) +
    1096      $           (-zdqsdif(ig)) *subtimestep) *z1(ig)
     1151     $           (-zdqsdif_tot(ig)) *subtimestep) *z1(ig)
    10971152                 zq1temp(ig)=zc(ig,1)
    10981153                 ENDIF !if h2o_ice_depth>0 and lag_layer
    … …  
    11071162     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
    11081163               !needs to change to the mean of eq
    1109                zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1)
     1164               zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1)
    11101165                !!!!! zdsqdif_ssi_frosst need to be changed to an
    11111166                !average
    … …  
    11151170     &                         /h2o_ice_depth(ig,1))
    11161171     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
    1117                zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1)
     1172               zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1)
    11181173               !needs to change to the mean of eq
    11191174              ENDIF
    … …  
    11341189                lh=(2834.3-0.28*(ztsrf(ig)-To)-
    11351190     &              0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3
    1136                 zdtsrf(ig,islope)=  zdqsdif(ig)*lh /pcapcal(ig,islope)
     1191                zdtsrf(ig,islope)=  zdqsdif_tot(ig)*lh
     1192     &                              /pcapcal(ig,islope)
    11371193              endif ! (latentheat_surfwater) then
    11381194
    … …  
    11461202     &                 + zdtsrf(ig,islope))*subtimestep
    11471203              zqsurf(ig)= zqsurf(ig)+(
    1148      &                       zdqsdif(ig))*subtimestep
     1204     &                       zdqsdif_tot(ig))*subtimestep
    11491205              if (zqsurf(ig)<0) then
    11501206                      zqsurf(ig)=0
    … …  
    11581214c             is still ice on the surface (oldschool)
    11591215               if(zqsurf(ig)
    1160      &           +zdqsdif(ig)*subtimestep
     1216     &           +zdqsdif_tot(ig)*subtimestep
    11611217     &           .gt.frost_albedo_threshold) then ! if there is still ice, T cannot exceed To
    11621218                 zdtsrf(ig,islope) = min(zdtsrf(ig,islope),
    … …  
    12971353     &                     "W.m-2",surf_h2o_lh(:,iflat))
    12981354         call write_output('zdq_subtimestep',
    1299      &                     'Actual flux zdqsdif*subtimestep',
    1300      &                     'kg.m-2',zdqsdif(:)*subtimestep)
     1355     &                     'Actual flux zdqsdif_surf*subtimestep',
     1356     &                     'kg.m-2',zdqsdif_surf(:)*subtimestep)
    13011357         call write_output('zdq_end',
    13021358     &                     'Flux after all contributions',
    1303      &                     'kg.m-2.s-1',zdqsdif(:))
     1359     &                     'kg.m-2.s-1',zdqsdif_surf(:))
    13041360C       Diagnostic output for HDO
    13051361!        if (hdo) then
Note: See TracChangeset for help on using the changeset viewer.