Ignore:
Timestamp:
Mar 13, 2020, 7:32:42 PM (5 years ago)
Author:
jnaar
Message:

[MARS GCM] Adding a new 2D variable "watercap", perennial ground water ice.
qsurf(ig,igcm_h2o_ice) can no longer be negative.
When watercaptag=true, ie when there is an infinite amount of water avalaible for sublimation,
and there is no more h2o frost, sublimation digs into watercap reservoir.
For now watercap can't grow, even on watercaptag, and has same albedo than qsurf(h2o_ice).
When using old start file, watercap is automatically initiated as 0.
Added watercap in newstart.F aswell.
JN

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

Legend:

Unmodified
Added
Removed

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

    r2257 r2260  
    840840                No=Nccnco2*tauscaling(ig)
    841841                Rn=-dlog(riceco2(ig,l))
    842                 n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
     842                n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
    843843                Qext1bins2(ig,l)=0.
    844844                do i = 1, nbinco2_cld
    845845                 n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
    846                  n_derf = erf((rb_cldco2(i+1)+Rn) *dev2)
     846                 n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
    847847                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf
    848848                 Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i)

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

    r2257 r2260  
    552552c Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance !
    553553             
    554               if ( (Ic_rice.ne.Ic_rice) ! will be true if it is Nan
    555      &               .or. Ic_rice .eq. 0.) then
     554              if (isnan(Ic_rice) .or. Ic_rice .eq. 0.) then
    556555                 Ic_rice=0.
    557556                 subpdtcloudco2(ig,l)=-sum_subpdt(ig,l)

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

    r2220 r2260  
    88                     day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf,co2ice, &
    99                     tauscaling,totcloudfrac,wstar,mem_Mccn_co2,mem_Nccn_co2,&
    10                      mem_Mh2o_co2)
     10                     mem_Mh2o_co2,watercap)
    1111
    1212  use tracer_mod, only: noms ! tracer names
    … …  
    6161  real,intent(out) :: mem_Nccn_co2(ngrid,nlay) ! Memory of CCN number of H2O and dust used by CO2
    6262  real,intent(out) :: mem_Mh2o_co2(ngrid,nlay) ! Memory of H2O mass integred into CO2 crystal
     63  real,intent(out) :: watercap(ngrid) ! h2o_ice_cover
    6364!======================================================================
    6465!  Local variables:
    … …  
    446447      write(*,*) 'phyetat0: loading surface tracer', &
    447448                           ' h2o_ice instead of h2o_vap'
     449      write(*,*) 'iq = ', iq
    448450    endif
    449451    call get_field(txt,qsurf(:,iq),found,indextime)
    … …  
    458460endif ! of if (nq.ge.1)
    459461
     462
     463! Surface water ice :
     464write(*,*) "qsurf(vap) : ", qsurf(:,2)
     465write(*,*) "qsurf(ice) : ", qsurf(:,3)
     466
     467
     468call get_field("watercap",watercap,found,indextime)
     469if (.not.found) then
     470  write(*,*) "phyetat0: Failed loading <watercap> : ", &
     471                       "<watercap> is set to zero"
     472  watercap(:)=0
     473else
     474  write(*,*) "phyetat0: Surface water ice <watercap> range:", &
     475             minval(watercap), maxval(watercap)
     476endif
     477
     478   write(*,*) 'Surface water ice cannot be allowed', &
     479                   ' to be negative if not on watercap !'
     480! tracer on surface
     481if (nq.ge.1) then
     482  do iq=1,nq
     483    txt=noms(iq)
     484!    if ((txt.eq."h2o_vap").or.(txt.eq."h2o_ice")) then
     485    if (txt.eq."h2o_ice") then
     486      do ig=1,ngrid
     487       if (qsurf(ig,iq).lt.0.0) then
     488          watercap(ig) = qsurf(ig,iq)
     489          qsurf(ig,iq) = 0.0
     490!        else if (qsurf(ig,iq)>0.5) then
     491!          watercap(ig) = watercap(ig) + qsurf(ig,iq) - 0.5
     492!          qsurf(ig,iq)=0.5
     493       end if
     494      end do
     495    endif
     496!    if (txt.eq."h2o_vap") then
     497!      do ig=1,ngrid
     498!        qsurf(ig,iq)=0.0
     499!        else if (qsurf(ig,iq)>0.5) then
     500!          watercap(ig) = watercap(ig) + qsurf(ig,iq) - 0.5
     501!          qsurf(ig,iq)=0.5
     502!      end do
     503!    end if
     504  enddo
     505endif
     506! tracer on surface
     507!if (nq.ge.1) then
     508!  do iq=1,nq
     509!    txt=noms(iq)
     510!    if (txt.eq."h2o_vap") then
     511      ! There is no surface tracer for h2o_vap;
     512      ! "h2o_ice" should be loaded instead
     513!      txt="h2o_ice"
     514!      write(*,*) 'phyetat0: loading surface tracer', &
     515!                           ' h2o_ice instead of h2o_vap'
     516!do ig=1,ngrid
     517! if (qsurf(ig,igcm_h2o_vap).lt.0.0) then
     518!    watercap(ig) = watercap(ig) + qsurf(ig,igcm_h2o_vap)
     519!    qsurf(ig,igcm_h2o_vap)=0.0
     520!  else if (qsurf(ig,iq)>0.5) then
     521!    watercap(ig) = watercap(ig) + qsurf(ig,iq) - 0.5
     522!    qsurf(ig,iq)=0.5
     523! end if
     524!end do
     525!    endif
     526!  enddo
     527!endif
     528
     529
     530
     531!do ig=1,ngrid
     532!  if (qsurf(ig,'h2o_ice').lt.0.0) then
     533!    watercap(ig) = watercap(ig) + qsurf(ig,'h2o_ice')
     534!    qsurf(ig,'h2o_ice')=0.0
     535!  else if (qsurf(ig,iq)>0.5) then
     536!    watercap(ig) = watercap(ig) + qsurf(ig,iq) - 0.5
     537!    qsurf(ig,iq)=0.5
     538!  end if
     539!end do
     540
    460541! Call to soil_settings, in order to read soil temperatures,
    461542! as well as thermal inertia and volumetric heat capacity

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

    r2220 r2260  
    150150                    phystep,time,tsurf,tsoil,co2ice,albedo,emis,q2,qsurf,&
    151151                    tauscaling,totcloudfrac,wstar, &
    152                     mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2)
     152                    mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2, watercap)
    153153  ! write time-dependent variable to restart file
    154154  use iostart, only : open_restartphy, close_restartphy, &
    … …  
    181181  real,intent(in) :: mem_Nccn_co2(ngrid,nlay) ! CCN number of H2O and dust used by CO2
    182182  real,intent(in) :: mem_Mh2o_co2(ngrid,nlay) ! H2O mass integred into CO2 crystal
     183  real,intent(in) :: watercap(ngrid)
    183184 
    184185  integer :: iq
    … …  
    198199  ! CO2 ice layer
    199200  call put_field("co2ice","CO2 ice cover",co2ice,time)
     201
     202  ! Water ice layer
     203  call put_field("watercap","H2O ice cover",watercap,time)
    200204 
    201205  ! Surface temperature

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

    r2252 r2260  
    4545     &                     tsurf, co2ice, emis,
    4646     &                     capcal, fluxgrd, qsurf,
    47      &                     hmons,summit,base
     47     &                     hmons,summit,base,watercap
    4848      use comsaison_h, only: dist_sol, declin, mu0, fract, local_time
    4949      use slope_mod, only: theta_sl, psi_sl
    … …  
    331331      real zdqmoldiff(ngrid,nlayer,nq)
    332332
     333      REAL dwatercap(ngrid), dwatercap_dif(ngrid)     ! (kg/m-2)
     334
    333335c     Local variable for local intermediate calcul:
    334336      REAL zflubid(ngrid)
    … …  
    512514     &         q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar,
    513515     &         mem_Mccn_co2,mem_Nccn_co2,
    514      &         mem_Mh2o_co2)
     516     &         mem_Mh2o_co2,watercap)
    515517
    516518         if (pday.ne.day_ini) then
    … …  
    539541      PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngrid,nq)
    540542      PRINT*,'check: co2 ',co2ice(1),co2ice(ngrid)
     543      PRINT*,'check: watercap ',watercap(1),watercap(ngrid)
    541544      !!!
    542545      day_ini = pday
    … …  
    675678      dsords(:,:)=0.
    676679      dsotop(:,:)=0.
     680      dwatercap(:)=0
    677681
    678682#ifdef DUSTSTORM
    … …  
    12411245     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
    12421246     &        zcondicea_co2microp,sensibFlux,
    1243      &        dustliftday,local_time)
     1247     &        dustliftday,local_time,watercap,dwatercap_dif)
    12441248
    12451249          DO ig=1,ngrid
    12461250             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
     1251             dwatercap(ig)=dwatercap(ig)+dwatercap_dif(ig)
    12471252          ENDDO
    12481253
    … …  
    20082013
    20092014c-----------------------------------------------------------------------
     2015c   J. Naar : Surface and sub-surface water ice
     2016c-----------------------------------------------------------------------
     2017c
     2018c
     2019c   Increment Watercap (surface h2o reservoirs):
     2020c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     2021
     2022      DO ig=1,ngrid
     2023         watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig)
     2024      ENDDO
     2025
     2026c     Test watercap
     2027c      DO ig=1,ngrid
     2028c         qsurf(iq,igcm_h2o_ice)=qsurf(iq,igcm_h2o_ice)+watercap(ig)
     2029c      ENDDO
     2030
     2031c-----------------------------------------------------------------------
    20102032c  13. Write output files
    20112033c  ----------------------
    … …  
    21852207     .                tsurf,tsoil,co2ice,albedo,emis,
    21862208     .                q2,qsurf,tauscaling,totcloudfrac,wstar,
    2187      .                mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2)
     2209     .                mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2,watercap)
    21882210         
    21892211         ENDIF
    … …  
    23682390        call wstats(ngrid,"co2ice","CO2 ice cover",
    23692391     &                "kg.m-2",2,co2ice)
     2392        call wstats(ngrid,"watercap","H2O ice cover",
     2393     &                "kg.m-2",2,watercap)
    23702394        call wstats(ngrid,"tauref","reference dod at 610 Pa","NU",
    23712395     &                2,tauref)
    … …  
    26882712         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
    26892713     &                                         ,"kg.m-2",2,co2ice)
     2714         call WRITEDIAGFI(ngrid,"watercap","Water ice thickness"
     2715     &                                         ,"kg.m-2",2,watercap)
    26902716
    26912717         call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",

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

    r2079 r2260  
    3131  REAL,SAVE,ALLOCATABLE :: fluxgrd(:) ! surface conduction flux (W.m-2)
    3232  REAL,ALLOCATABLE,SAVE :: qsurf(:,:) ! tracer on surface (e.g. kg.m-2)
     33  REAL,SAVE,ALLOCATABLE :: watercap(:) ! Surface water ice (kg.m-2)
    3334
    3435contains
    … …  
    5354    allocate(tsurf(ngrid))
    5455    allocate(co2ice(ngrid))
     56    allocate(watercap(ngrid))
    5557    allocate(emis(ngrid))
    5658    allocate(capcal(ngrid))
    … …  
    8082    if (allocated(tsurf))       deallocate(tsurf)
    8183    if (allocated(co2ice))      deallocate(co2ice)
     84    if (allocated(watercap))    deallocate(watercap)
    8285    if (allocated(emis))        deallocate(emis)
    8386    if (allocated(capcal))      deallocate(capcal)

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

    r2218 r2260  
    1313     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
    1414     $                hfmax,pcondicea_co2microp,sensibFlux,
    15      $                dustliftday,local_time)
     15     $                dustliftday,local_time,watercap, dwatercap_dif)
    1616
    1717      use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
    … …  
    142142      REAL kmixmin
    143143
     144c    Argument added for surface water ice budget:
     145      REAL,INTENT(IN) :: watercap(ngrid)
     146      REAL,INTENT(OUT) :: dwatercap_dif(ngrid)
     147
    144148c    Mass-variation scheme :
    145149c    ~~~~~~~
    … …  
    213217      pdqdif(1:ngrid,1:nlay,1:nq)=0
    214218      pdqsdif(1:ngrid,1:nq)=0
     219      dwatercap_dif(1:ngrid)=0
    215220
    216221c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
    … …  
    838843
    839844            DO ig=1,ngrid
    840               if(.not.watercaptag(ig)) then
    841                 if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
    842      &             .gt.pqsurf(ig,igcm_h2o_ice)) then
    843 c                 write(*,*)'on sublime plus que qsurf!'
    844                   pdqsdif(ig,igcm_h2o_ice)=
     845              if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
     846     &             .gt.(pqsurf(ig,igcm_h2o_ice))) then
     847c              write(*,*)'subliming more than available frost:  qsurf!'
     848                if(watercaptag(ig)) then
     849c              dwatercap_dif same sign as pdqsdif (pos. toward ground)
     850                   dwatercap_dif(ig) = -pdqsdif(ig,igcm_h2o_ice)
     851     &                         - pqsurf(ig,igcm_h2o_ice)/ptimestep
     852                   pdqsdif(ig,igcm_h2o_ice)=
    845853     &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
     854c                  write(*,*) "subliming more than available frost:  qsurf!"
     855                else  !  (case not.watercaptag(ig))
     856                   pdqsdif(ig,igcm_h2o_ice)=
     857     &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
     858                   dwatercap_dif(ig) = 0.0
     859
     860                endif   
    846861c                 write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
    847862                  z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
    848863                  zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
    849      $            zb(ig,2)*zc(ig,2) +
    850      $            (-pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig)
     864     $            zb(ig,2)*zc(ig,2) + (dwatercap_dif(ig)-
     865c     $            zb(ig,2)*zc(ig,2) + (-
     866     $            pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig)
    851867                  zq1temp(ig)=zc(ig,1)
    852                 endif   
    853868              endif ! if (.not.watercaptag(ig))
    854869c             Starting upward calculations for water :
    … …  
    871886                lh=(2834.3-0.28*(tsrf_lh(ig)-To)-
    872887     &              0.004*(tsrf_lh(ig)-To)*(tsrf_lh(ig)-To))*1.e+3
    873                 pdtsrf(ig)= pdtsrf(ig) + pdqsdif(ig,igcm_h2o_ice)*lh
    874      &                   /pcapcal(ig)
     888                pdtsrf(ig)= pdtsrf(ig) + (dwatercap_dif(ig)+
     889     &               pdqsdif(ig,igcm_h2o_ice))*lh /pcapcal(ig)
    875890              endif
    876891
    877892               if(pqsurf(ig,igcm_h2o_ice)
    878      &           +pdqsdif(ig,igcm_h2o_ice)*ptimestep
     893     &           +(dwatercap_dif(ig)+pdqsdif(ig,igcm_h2o_ice))*ptimestep
    879894     &           .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To
    880895     &           pdtsrf(ig) = min(pdtsrf(ig),(To-ptsrf(ig))/ptimestep) ! ice melt case
Note: See TracChangeset for help on using the changeset viewer.